Warning, /GeneratorInterface/PartonShowerVeto/src/ME2pythia.f_old is written in an unsupported language. File is not indexed.
0001 C******************************************************
0002 C* MadEvent - Pythia interface. *
0003 C* Version 4.2, 4 March 2007 *
0004 C* *
0005 C* dkcira 2008.02.05 *
0006 C* *
0007 C* - Improvement of matching routines *
0008 C* *
0009 C* Version 4.1 *
0010 C* *
0011 C* - Possibility to use several files *
0012 C* *
0013 C* Version 4.0 *
0014 C* *
0015 C* - Routines for matching of ME and PS *
0016 C* *
0017 C* Version 3.8 *
0018 C* *
0019 C* - Give the event number in the event file in the *
0020 C* new variable IEVNT in UPPRIV *
0021 C* *
0022 C* Version 3.7 *
0023 C* *
0024 C* - Set mass of massless outgoing particles to *
0025 C* Pythia mass (PMAS(I,1)) *
0026 C* *
0027 C* Version 3.6 *
0028 C* *
0029 C* - Removed the 1st # from the event file header *
0030 C* *
0031 C* Version 3.5 *
0032 C* *
0033 C* - Reads according to the new LH event file format *
0034 C* - Now only LNHIN, LNHOUT and MSCAL in UPPRIV *
0035 C* *
0036 C* Version 3.4 *
0037 C* *
0038 C* - Reads particle masses from event file *
0039 C* *
0040 C* Version 3.3 *
0041 C* *
0042 C* - Added option MSCAL in common block UPPRIV to *
0043 C* choose between fix (0) or event-based (1) *
0044 C* scale for Pythia parton showering (SCALUP). *
0045 C* - Fixed bug in reading the SLHA file *
0046 C* *
0047 C* Version 3.2 *
0048 C* *
0049 C* - Reading the SLHA format param_card from the *
0050 C* banner *
0051 C* - Added support for lpp1/lpp2 = 2 or 3 *
0052 C* - Removed again support for different MadEvent *
0053 C* processes in different files (no longer *
0054 C* necessary with new multiple processes support *
0055 C* in MadGraph/MadEvent *
0056 C* *
0057 C* Version 3.1 *
0058 C* - Added support for different MadEvent processes *
0059 C* in different files *
0060 C* - Fixed bug in e+e- collisions *
0061 C* *
0062 C* Written by J.Alwall, alwall@fyma.ucl.ac.be *
0063 C* Earlier versions by S.Mrenna, M.Kirsanov *
0064 C* *
0065 C******************************************************
0066 C* *
0067 C* Instructions: *
0068 C* Please use the common block UPPRIV: *
0069 C* - The logical unit LNHIN must be an opened *
0070 C* MadEvent event file *
0071 C* - The output unit LNHOUT is by default 6 (std out) *
0072 C* - Set MSCAL to 1 if a dynamical scale is desired *
0073 C* for parton showers rather than the one given as *
0074 C* factorization scale by MadEvent (otherwise 0) *
0075 C* - IEVNT gives the number of the event in the event *
0076 C* file *
0077 C* - ICKKW is set automatically depending on whether *
0078 C* the events generated are matched or not *
0079 C* *
0080 C******************************************************
0081
0082 C*********************************************************************
0083 C...UPINIT
0084 C...Routine called by PYINIT to set up user-defined processes.
0085 C*********************************************************************
0086 SUBROUTINE MGINIT(npara,param,value)
0087
0088 IMPLICIT NONE
0089 c
0090 c arguments
0091 c
0092 integer npara
0093 character*20 param(*),value(*)
0094
0095 C CHARACTER*132 CHAR_READ
0096
0097 C...Pythia parameters.
0098 INTEGER MSTP,MSTI,MRPY
0099 DOUBLE PRECISION PARP,PARI,RRPY
0100 COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0101 COMMON/PYDATR/MRPY(6),RRPY(100)
0102
0103 C...User process initialization commonblock.
0104 INTEGER MAXPUP
0105 PARAMETER (MAXPUP=100)
0106 INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
0107 DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
0108 COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
0109 & IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP),
0110 & LPRUP(MAXPUP)
0111
0112 C...Extra commonblock to transfer run info.
0113 INTEGER LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
0114 COMMON/UPPRIV/LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
0115 DATA LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE/77,6,1,0,0,1/
0116 SAVE/UPPRIV/
0117
0118 C...Inputs for the matching algorithm
0119 double precision etcjet,rclmax,etaclmax,qcut,clfact
0120 integer maxjets,minjets,iexcfile,ktsche,nexcres,excres(30)
0121 common/MEMAIN/etcjet,rclmax,etaclmax,qcut,clfact,
0122 $ maxjets,minjets,iexcfile,ktsche,nexcres,excres
0123
0124 C
0125 C...Set kt clustering scheme (if not already set)
0126 C
0127 IF(ABS(IDBMUP(1)).EQ.11.AND.ABS(IDBMUP(2)).EQ.11.AND.
0128 $ IDBMUP(1).EQ.-IDBMUP(2).AND.ktsche.EQ.0)THEN
0129 ktsche=1
0130 ELSE IF(ktsche.EQ.0) THEN
0131 ktsche=4313
0132 ENDIF
0133
0134 IF(ickkw.gt.0) CALL set_matching(npara,param,value)
0135
0136 C...For photon initial states from protons: Set proton not to break up
0137 CALL PYGIVE('MSTP(98)=1')
0138
0139 RETURN
0140
0141 END
0142
0143 C*********************************************************************
0144 C...UPEVNT
0145 C...Routine called by PYEVNT or PYEVNW to get user process event
0146 C*********************************************************************
0147 SUBROUTINE MGEVNT
0148
0149 IMPLICIT NONE
0150
0151 C...Pythia parameters.
0152 INTEGER MSTP,MSTI
0153 DOUBLE PRECISION PARP,PARI
0154 COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0155
0156 C...User process initialization commonblock.
0157 INTEGER MAXPUP
0158 PARAMETER (MAXPUP=100)
0159 INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
0160 DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
0161 COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
0162 & IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP),
0163 & LPRUP(MAXPUP)
0164
0165 C...User process event common block.
0166 INTEGER MAXNUP
0167 PARAMETER (MAXNUP=500)
0168 INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
0169 DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
0170 COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP),
0171 & ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP),
0172 & VTIMUP(MAXNUP),SPINUP(MAXNUP)
0173
0174 C...Pythia common blocks
0175 INTEGER PYCOMP,KCHG,MINT,NPART,NPARTD,IPART,MAXNUR
0176 DOUBLE PRECISION PMAS,PARF,VCKM,VINT,PTPART
0177
0178 C...Particle properties + some flavour parameters.
0179 COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0180 COMMON/PYINT1/MINT(400),VINT(400)
0181 PARAMETER (MAXNUR=1000)
0182 COMMON/PYPART/NPART,NPARTD,IPART(MAXNUR),PTPART(MAXNUR)
0183
0184 C...Extra commonblock to transfer run info.
0185 INTEGER LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
0186 COMMON/UPPRIV/LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
0187
0188 C...Inputs for the matching algorithm
0189 double precision etcjet,rclmax,etaclmax,qcut,clfact
0190 integer maxjets,minjets,iexcfile,ktsche,nexcres,excres(30)
0191 common/MEMAIN/etcjet,rclmax,etaclmax,qcut,clfact,
0192 $ maxjets,minjets,iexcfile,ktsche,nexcres,excres
0193
0194 C...Commonblock to transfer event-by-event matching info
0195 INTEGER NLJETS,IEXC,Ifile
0196 DOUBLE PRECISION PTCLUS
0197 COMMON/MEMAEV/PTCLUS(20),NLJETS,IEXC,Ifile
0198
0199 C...Local variables
0200 INTEGER I,J,IBEG,NEX,KP(MAXNUP),MOTH,NUPREAD!,IREM
0201 DOUBLE PRECISION PSUM,ESUM,PM1,PM2,A1,A2,A3,A4,A5
0202 DOUBLE PRECISION SCALLOW(MAXNUP),PNONJ(4),PMNONJ!,PT2JETS
0203
0204 NUPREAD=NUP
0205
0206 C...Read NUP subsequent lines with information on each particle.
0207 115 ESUM=0d0
0208 PSUM=0d0
0209 NEX=2
0210 NUP=1
0211
0212 DO 120 I=1,NUPREAD
0213 C...Reset resonance momentum to prepare for mass shifts
0214 IF(ISTUP(NUP).EQ.2) PUP(3,NUP)=0
0215 IF(ISTUP(NUP).EQ.1)THEN
0216 NEX=NEX+1
0217 IF(PUP(5,NUP).EQ.0D0.AND.IABS(IDUP(NUP)).GT.3
0218 $ .AND.IDUP(NUP).NE.21) THEN
0219 C...Set massless particle masses to Pythia default. Adjust z-momentum.
0220 PUP(5,NUP)=PMAS(IABS(PYCOMP(IDUP(NUP))),1)
0221 PUP(3,NUP)=SIGN(SQRT(MAX(0d0,PUP(4,NUP)**2-PUP(5,NUP)**2-
0222 $ PUP(1,NUP)**2-PUP(2,NUP)**2)),PUP(3,NUP))
0223 ENDIF
0224 PSUM=PSUM+PUP(3,NUP)
0225 C...Set mother resonance momenta
0226 MOTH=MOTHUP(1,NUP)
0227 DO WHILE (MOTH.GT.2)
0228 PUP(3,MOTH)=PUP(3,MOTH)+PUP(3,NUP)
0229 MOTH=MOTHUP(1,MOTH)
0230 ENDDO
0231 ENDIF
0232 NUP=NUP+1
0233 120 CONTINUE
0234 NUP=NUP-1
0235
0236 C..Adjust mass of resonances
0237 DO I=1,NUP
0238 IF(ISTUP(I).EQ.2)THEN
0239 PUP(5,I)=SQRT(PUP(4,I)**2-PUP(1,I)**2-PUP(2,I)**2-
0240 $ PUP(3,I)**2)
0241 ENDIF
0242 ENDDO
0243
0244 C...Adjust energy and momentum of incoming particles
0245 C...In massive case need to solve quadratic equation
0246 c PM1=PUP(5,1)**2
0247 c PM2=PUP(5,2)**2
0248 c A1=4d0*(ESUM**2-PSUM**2)
0249 c A2=ESUM**2-PSUM**2+PM2-PM1
0250 c A3=2d0*PSUM*A2
0251 c A4=A3/A1
0252 c A5=(A2**2-4d0*ESUM**2*PM2)/A1
0253 c
0254 c PUP(3,2)=A4+SIGN(SQRT(A4**2+A5),PUP(3,2))
0255 c PUP(3,1)=PSUM-PUP(3,2)
0256 c PUP(4,1)=SQRT(PUP(3,1)**2+PM1)
0257 c PUP(4,2)=SQRT(PUP(3,2)**2+PM2)
0258
0259 ESUM=PUP(4,1)+PUP(4,2)
0260
0261 C...Assuming massless incoming particles - otherwise Pythia adjusts
0262 C...the momenta to make them massless
0263 IF(IDBMUP(1).GT.100.AND.IDBMUP(2).GT.100)THEN
0264 DO I=1,2
0265 PUP(3,I)=0.5d0*(PSUM+SIGN(ESUM,PUP(3,I)))
0266 PUP(5,I)=0d0
0267 ENDDO
0268 PUP(4,1)=ABS(PUP(3,1))
0269 PUP(4,2)=ESUM-PUP(4,1)
0270 ENDIF
0271
0272 C...If you want to use some other scale for parton showering then the
0273 C...factorisation scale given by MadEvent, please implement the function PYMASC
0274 C...(example function included below)
0275
0276 IF(ickkw.eq.0.AND.MSCAL.GT.0) CALL PYMASC(SCALUP)
0277 IF(MINT(35).eq.3.AND.ickkw.EQ.1) SCALUP=SQRT(PARP(67))*SCALUP
0278
0279 IF(ickkw.gt.0) THEN
0280 c
0281 c Set up number of jets
0282 c
0283 NLJETS=0
0284 NPART=0
0285 do i=3,NUP
0286 if(ISTUP(i).ne.1) cycle
0287 NPART=NPART+1
0288 IPART(NPART)=i
0289 if(iabs(IDUP(i)).ge.6.and.IDUP(i).ne.21) cycle
0290 if(MOTHUP(1,i).gt.2) cycle
0291 NLJETS=NLJETS+1
0292 PTCLUS(NLJETS)=PTPART(NPART)
0293 enddo
0294 CALL ALPSOR(PTCLUS,nljets,KP,1)
0295
0296 C IDPRUP=LPRUP(NLJETS-MINJETS+1)
0297
0298 C IF(ickkw.eq.1) THEN
0299 Cc ... and decide whether exclusive or inclusive
0300 C if(IEXCFILE.eq.0.and.NLJETS.eq.MAXJETS)then
0301 C IEXC=0
0302 C else
0303 C IEXC=1
0304 C endif
0305 C ENDIF
0306 ENDIF
0307
0308 RETURN
0309 END
0310
0311 C*********************************************************************
0312 C...UPVETO
0313 C...Subroutine to implement the MLM jet matching criterion
0314 C*********************************************************************
0315 SUBROUTINE MGVETO(IPVETO)
0316
0317 IMPLICIT NONE
0318
0319 C...Pythia common blocks
0320 INTEGER MINT
0321 DOUBLE PRECISION VINT
0322 COMMON/PYINT1/MINT(400),VINT(400)
0323
0324 C...GUP Event common block
0325 INTEGER MAXNUP
0326 PARAMETER (MAXNUP=500)
0327 INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
0328 DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
0329 COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
0330 & IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP),
0331 & ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP),
0332 & SPINUP(MAXNUP)
0333
0334 C...User process initialization commonblock.
0335 INTEGER MAXPUP
0336 PARAMETER (MAXPUP=100)
0337 INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
0338 DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
0339 COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
0340 & IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP),
0341 & LPRUP(MAXPUP)
0342
0343 C...HEPEVT commonblock.
0344 INTEGER NMXHEP,NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
0345 PARAMETER (NMXHEP=4000)
0346 COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
0347 &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
0348 DOUBLE PRECISION PHEP,VHEP
0349 SAVE /HEPEVT/
0350 INTEGER IPVETO
0351
0352 C...GETJET commonblocks
0353 INTEGER MNCY,MNCPHI,NCY,NCPHI,NJMAX,JETNO,NCJET
0354 DOUBLE PRECISION YCMIN,YCMAX,DELY,DELPHI,ET,STHCAL,CTHCAL,CPHCAL,
0355 & SPHCAL,PCJET,ETJET
0356 PARAMETER (MNCY=200)
0357 PARAMETER (MNCPHI=200)
0358 COMMON/CALORM/DELY,DELPHI,ET(MNCY,MNCPHI),
0359 $CTHCAL(MNCY),STHCAL(MNCY),CPHCAL(MNCPHI),SPHCAL(MNCPHI),
0360 $YCMIN,YCMAX,NCY,NCPHI
0361 DATA NCY,NCPHI/50,60/
0362 PARAMETER (NJMAX=500)
0363 COMMON/GETCOMM/PCJET(4,NJMAX),ETJET(NJMAX),JETNO(MNCY,MNCPHI),
0364 $NCJET
0365 DOUBLE PRECISION PI
0366 PARAMETER (PI=3.141593D0)
0367 C
0368 DOUBLE PRECISION PSERAP
0369 INTEGER K(NJMAX),KP(NJMAX),kpj(njmax)
0370
0371 C...Variables for the kT-clustering
0372 INTEGER NMAX,NN,NJET,NSUB,JET,NJETM,IHARD,IP1,IP2
0373 DOUBLE PRECISION PP,PJET
0374 DOUBLE PRECISION ECUT,Y,YCUT,RAD
0375 PARAMETER (NMAX=512)
0376 DIMENSION JET(NMAX),Y(NMAX),PP(4,NMAX),PJET(4,NMAX),
0377 $ PJETM(4,NMAX)
0378 INTEGER NNM
0379 DOUBLE PRECISION YM(NMAX),PPM(4,NMAX),PJETM
0380
0381 C...kt clustering common block
0382 INTEGER NMAXKT,NUM,HIST
0383 PARAMETER (NMAXKT=512)
0384 DOUBLE PRECISION PPP,KT,ETOT,RSQ,KTP,KTS,KTLAST
0385 COMMON /KTCOMM/ETOT,RSQ,PPP(9,NMAXKT),KTP(NMAXKT,NMAXKT),
0386 $ KTS(NMAXKT),KT(NMAXKT),KTLAST(NMAXKT),HIST(NMAXKT),NUM
0387
0388 C...Extra commonblock to transfer run info.
0389 INTEGER LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
0390 COMMON/UPPRIV/LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
0391
0392 C...Inputs for the matching algorithm
0393 C clfact determines how jet-to parton matching is done
0394 C kt-jets: default=1
0395 C clfact >= 0: Max mult. if within clfact*max(qcut,Q(partNmax)) from jet, others within clfact*qcut
0396 C clfact < 0: Max mult. if within |clfact|*Q(jetNmax) from jet, other within |clfact|*qcut
0397 C cone-jets: default=1.5
0398 C Matching if within clfact*RCLMAX
0399
0400 double precision etcjet,rclmax,etaclmax,qcut,clfact
0401 integer maxjets,minjets,iexcfile,ktsche,nexcres,excres(30)
0402 common/MEMAIN/etcjet,rclmax,etaclmax,qcut,clfact,
0403 $ maxjets,minjets,iexcfile,ktsche,nexcres,excres
0404
0405 C...Commonblock to transfer event-by-event matching info
0406 INTEGER NLJETS,IEXC,Ifile
0407 DOUBLE PRECISION PTCLUS
0408 COMMON/MEMAEV/PTCLUS(20),NLJETS,IEXC,Ifile
0409
0410 INTEGER nvarev,nvar2
0411 PARAMETER (nvarev=57,nvar2=6)
0412
0413 REAL*4 varev(nvarev)
0414 COMMON/HISTDAT/varev
0415
0416
0417 C local variables
0418 integer i,j,ihep,nmatch,jrmin,KPT(MAXNUP),nres
0419 double precision etajet,phijet,delr,dphi,delrmin,ptjet
0420 double precision p(4,10),pt(10),eta(10),phi(10)
0421 integer idbg
0422 data idbg/0/
0423 REAL*4 var2(nvar2)
0424 c
0425 double precision tiny
0426 parameter (tiny=1d-3)
0427 integer icount
0428 data icount/0/
0429 INTEGER ISTOLD(NMXHEP),IST,IDA
0430 c LOGICAL FOUND
0431
0432 c if(NLJETS.GT.0)then
0433 c idbg=1
0434 c else
0435 c idbg=0
0436 c endif
0437
0438 IPVETO=0
0439
0440 IF(ICKKW.NE.1) RETURN
0441
0442 IF(NLJETS.LT.MINJETS.OR.NLJETS.GT.MAXJETS)THEN
0443 if(idbg.eq.1)
0444 $ WRITE(LNHOUT,*) 'Failed due to NLJETS ',NLJETS,' < ',MINJETS,
0445 $ ' or > ',MAXJETS
0446 GOTO 999
0447 ENDIF
0448
0449 C Throw event if it contains an excluded resonance
0450 NRES=0
0451 DO I=1,NUP
0452 IF(ISTUP(I).EQ.2)THEN
0453 DO J=1,nexcres
0454 IF(IDUP(I).EQ.EXCRES(J)) NRES=NRES+1
0455 ENDDO
0456 ENDIF
0457 ENDDO
0458 IF(NRES.GT.0)THEN
0459 if(idbg.eq.1)
0460 $ PRINT *,'Event',IEVNT,' thrown because of ',NRES,
0461 $ ' excluded resonance(s)'
0462 c CALL PYLIST(7)
0463 GOTO 999
0464 ENDIF
0465
0466
0467
0468 C Set up vetoed mothers
0469 c DO I=1,MAXNUP
0470 c INORAD(I)=0
0471 c ENDDO
0472 c DO IHEP=1,NUP-2
0473 c if(ISTHEP(ihep).gt.1.and.iabs(IDHEP(ihep)).gt.8) then
0474 c if(iabs(IDHEP(ihep)).gt.5.and.IDHEP(ihep).ne.21) then
0475 c INORAD(ihep)=1
0476 c endif
0477 c ENDDO
0478
0479 C Set status for non-clustering partons to 2
0480 DO ihep=1,NHEP
0481 c ISTORG(ihep)=ISTHEP(ihep)
0482 IF(ISTHEP(ihep).EQ.1.AND.iabs(IDHEP(ihep)).GE.6.AND.
0483 $ IDHEP(ihep).NE.21) ISTHEP(ihep)=2
0484 IF(ISTHEP(ihep).EQ.1.AND.(iabs(IDHEP(ihep)).lt.6.or.
0485 $ IDHEP(ihep).eq.21).AND.JMOHEP(1,ihep).GT.0) then
0486 IDA=ihep
0487 DO WHILE(JMOHEP(1,IDA).GT.0)
0488 c Trace mothers, if daughter particle in hard event,
0489 c must be decay - remove
0490 IF(iabs(IDHEP(JMOHEP(1,IDA))).GE.6.AND.
0491 $ IDHEP(JMOHEP(1,IDA)).NE.21.AND.
0492 $ IDA.LE.NUP+4) GOTO 5
0493 IDA=JMOHEP(1,IDA)
0494 ENDDO
0495 cycle
0496 5 ISTHEP(ihep)=2
0497 ENDIF
0498 ENDDO
0499
0500 C Prepare histogram filling
0501 DO I=1,4
0502 var2(1+I)=-1
0503 varev(46+I)=-1
0504 varev(50+I)=-1
0505 ENDDO
0506
0507 C CHECK FOR EVENT ERROR OR ZERO WGT
0508 I=0
0509 C
0510 c reconstruct parton-level event
0511
0512 if(idbg.eq.1) then
0513 write(LNHOUT,*) ' '
0514 write(LNHOUT,*) 'new event '
0515 c CALL PYLIST(1)
0516 CALL PYLIST(7)
0517 CALL PYLIST(5)
0518 write(LNHOUT,*) 'PARTONS'
0519 endif
0520 i=0
0521 do ihep=3,nup
0522 if(ISTUP(ihep).ne.1.or.
0523 $ (iabs(IDUP(ihep)).ge.6.and.IDUP(ihep).ne.21)) cycle
0524 if(MOTHUP(1,ihep).gt.2) cycle
0525 i=i+1
0526 do j=1,4
0527 p(j,i)=pup(j,ihep)
0528 enddo
0529 pt(i)=sqrt(p(1,i)**2+p(2,i)**2)
0530 if(i.LE.4) varev(50+i)=pt(i)
0531 eta(i)=-log(tan(0.5d0*atan2(pt(i)+tiny,p(3,i))))
0532 phi(i)=atan2(p(2,i),p(1,i))
0533 if(idbg.eq.1) then
0534 write(LNHOUT,*) pt(i),eta(i),phi(i)
0535 endif
0536 enddo
0537 if(i.ne.NLJETS)then
0538 print *,'Error in UPVETO: Wrong number of jets found ',i,NLJETS
0539 CALL PYLIST(7)
0540 CALL PYLIST(2)
0541 stop
0542 endif
0543
0544 C Bubble-sort PTs in descending order
0545 DO I=1,3
0546 DO J=4,I+1,-1
0547 IF(varev(50+J).GT.varev(50+I))THEN
0548 PTJET=varev(50+J)
0549 varev(50+J)=varev(50+I)
0550 varev(50+I)=PTJET
0551 ENDIF
0552 ENDDO
0553 ENDDO
0554
0555 if(idbg.eq.1) then
0556 do i=1,nhep
0557 write(LNHOUT,111) i,isthep(i),idhep(i),jmohep(1,i),jmohep(2,i)
0558 $ ,phep(1,i),phep(2,i),phep(3,i)
0559 enddo
0560 111 format(5(i4,1x),3(f12.5,1x))
0561 endif
0562
0563 IF(qcut.le.0d0)then
0564
0565 IF(clfact.EQ.0d0) clfact=1.5d0
0566
0567 c CALL PYLIST(7)
0568 c CALL PYLIST(2)
0569 c CALL PYLIST(5)
0570 c Start from the partonic system
0571 IF(NLJETS.GT.0) CALL ALPSOR(pt,nljets,KP,2)
0572 c reconstruct showered jets
0573 c
0574 YCMAX=ETACLMAX+RCLMAX
0575 YCMIN=-YCMAX
0576 CALL CALINIM
0577 CALL CALDELM(1,1)
0578 CALL GETJETM(RCLMAX,ETCJET,ETACLMAX)
0579 c analyse only events with at least nljets-reconstructed jets
0580 IF(NCJET.GT.0) CALL ALPSOR(ETJET,NCJET,K,2)
0581 if(idbg.eq.1) then
0582 write(LNHOUT,*) 'JETS'
0583 do i=1,ncjet
0584 j=k(ncjet+1-i)
0585 ETAJET=PSERAP(PCJET(1,j))
0586 PHIJET=ATAN2(PCJET(2,j),PCJET(1,j))
0587 write(LNHOUT,*) etjet(j),etajet,phijet
0588 enddo
0589 endif
0590 IF(NCJET.LT.NLJETS) THEN
0591 if(idbg.eq.1)
0592 $ WRITE(LNHOUT,*) 'Failed due to NCJET ',NCJET,' < ',NLJETS
0593 GOTO 999
0594 endif
0595 c associate partons and jets, using min(delr) as criterion
0596 NMATCH=0
0597 DO I=1,NCJET
0598 KPJ(I)=0
0599 ENDDO
0600 DO I=1,NLJETS
0601 DELRMIN=1D5
0602 DO 110 J=1,NCJET
0603 IF(KPJ(J).NE.0) GO TO 110
0604 ETAJET=PSERAP(PCJET(1,J))
0605 PHIJET=ATAN2(PCJET(2,J),PCJET(1,J))
0606 DPHI=ABS(PHI(KP(NLJETS-I+1))-PHIJET)
0607 IF(DPHI.GT.PI) DPHI=2.*PI-DPHI
0608 DELR=SQRT((ETA(KP(NLJETS-I+1))-ETAJET)**2+(DPHI)**2)
0609 IF(DELR.LT.DELRMIN) THEN
0610 DELRMIN=DELR
0611 JRMIN=J
0612 ENDIF
0613 110 CONTINUE
0614 IF(DELRMIN.LT.clfact*RCLMAX) THEN
0615 NMATCH=NMATCH+1
0616 KPJ(JRMIN)=I
0617 ENDIF
0618 C WRITE(*,*) 'PARTON-JET',I,' best match:',k(ncjet+1-jrmin)
0619 c $ ,delrmin
0620 ENDDO
0621 IF(NMATCH.LT.NLJETS) THEN
0622 if(idbg.eq.1)
0623 $ WRITE(LNHOUT,*) 'Failed due to NMATCH ',NMATCH,' < ',NLJETS
0624 GOTO 999
0625 endif
0626 C REJECT EVENTS WITH LARGER JET MULTIPLICITY FROM EXCLUSIVE SAMPLE
0627 IF(NCJET.GT.NLJETS.AND.IEXC.EQ.1) THEN
0628 if(idbg.eq.1)
0629 $ WRITE(LNHOUT,*) 'Failed due to NCJET ',NCJET,' > ',NLJETS
0630 GOTO 999
0631 endif
0632 C VETO EVENTS WHERE MATCHED JETS ARE SOFTER THAN NON-MATCHED ONES
0633 IF(IEXC.NE.1) THEN
0634 J=NCJET
0635 DO I=1,NLJETS
0636 IF(KPJ(K(J)).EQ.0) GOTO 999
0637 J=J-1
0638 ENDDO
0639 ENDIF
0640
0641 else ! qcut.gt.0
0642
0643 if(MINT(35).eq.3) then
0644 C The pt-ordered showers have been used - use "shower emission pt method"
0645 C Veto events where first shower emission has kt > YCUT
0646
0647 IF(NLJETS.EQ.0)THEN
0648 VINT(358)=0
0649 ENDIF
0650
0651 IF(idbg.eq.1) THEN
0652 PRINT *,'Using shower emission pt method'
0653 PRINT *,'qcut, ptclus(1), vint(357),vint(358): ',
0654 $ qcut,ptclus(1),vint(357),vint(358)
0655 ENDIF
0656 YCUT=qcut**2
0657
0658 IF(NLJETS.GT.0.AND.PTCLUS(1)**2.LT.YCUT) THEN
0659 if(idbg.eq.1)
0660 $ WRITE(LNHOUT,*) 'Failed due to KT ',
0661 $ PTCLUS(1),' < ',SQRT(YCUT)
0662 GOTO 999
0663 ENDIF
0664
0665 c PRINT *,'Y,VINT:',SQRT(Y(NLJETS+1)),SQRT(VINT(390))
0666
0667 IF(IEXC.EQ.1.AND.MAX(VINT(357),VINT(358)).GT.SQRT(YCUT))THEN
0668 if(idbg.eq.1)
0669 $ WRITE(LNHOUT,*),
0670 $ 'Failed due to ',max(VINT(357),VINT(358)),' > ',SQRT(YCUT)
0671 GOTO 999
0672 ENDIF
0673 c PRINT *,NLJETS,IEXC,SQRT(VINT(390)),PTCLUS(1),SQRT(YCUT)
0674 c Highest multiplicity case
0675 IF(IEXC.EQ.0.AND.NLJETS.GT.0.AND.
0676 $ MAX(VINT(357),VINT(358)).GT.PTCLUS(1))THEN
0677 c $ VINT(390).GT.PTCLUS(1)**2)THEN
0678 if(idbg.eq.1)
0679 $ WRITE(LNHOUT,*),
0680 $ 'Failed due to ',max(VINT(357),VINT(358)),' > ',PTCLUS(1)
0681 GOTO 999
0682 ENDIF
0683 c
0684 else ! not false ! not pt-ordered showers
0685
0686 IF(clfact.EQ.0d0) clfact=1d0
0687
0688 C---FIND FINAL STATE COLOURED PARTICLES
0689 NN=0
0690 DO IHEP=1,NHEP
0691 IF (ISTHEP(IHEP).EQ.1
0692 & .AND.(ABS(IDHEP(IHEP)).LT.6.OR.IDHEP(IHEP).EQ.21)) THEN
0693 PTJET=sqrt(PHEP(1,IHEP)**2+PHEP(2,IHEP)**2)
0694 ETAJET=ABS(LOG(MIN((SQRT(PTJET**2+PHEP(3,IHEP)**2)+
0695 $ ABS(PHEP(3,IHEP)))/PTJET,1d5)))
0696 IF(ETAJET.GT.etaclmax) cycle
0697 NN=NN+1
0698 IF (NN.GT.NMAX) then
0699 CALL PYLIST(2)
0700 PRINT *, 'Too many particles: ', NN
0701 NN=NN-1
0702 GOTO 10
0703 endif
0704 DO I=1,4
0705 PP(I,NN)=PHEP(I,IHEP)
0706 ENDDO
0707 ELSE if(idbg.eq.1)THEN
0708 PRINT *,'Skipping particle ',IHEP,ISTHEP(IHEP),IDHEP(IHEP)
0709 ENDIF
0710 ENDDO
0711
0712 C...Cluster event to find values of Y including jet matching but not veto of too many jets
0713 C...Only used to fill the beforeveto Root tree
0714 10 ECUT=1
0715 IF (NN.GT.1) then
0716 CALL KTCLUS(KTSCHE,PP,NN,ECUT,Y,*999)
0717 if(idbg.eq.1)
0718 $ WRITE(LNHOUT,*) 'Clustering values:',
0719 $ (SQRT(Y(i)),i=1,MIN(NN,3))
0720
0721 C Print out values in the case where all jets are matched at the
0722 C value of the NLJETS:th clustering
0723 var2(1)=NLJETS
0724 var2(6)= Ifile
0725
0726 if(NLJETS.GT.MINJETS)then
0727 YCUT=Y(NLJETS)
0728 CALL KTRECO(MOD(KTSCHE,10),PP,NN,ECUT,YCUT,YCUT,PJET,JET,
0729 $ NCJET,NSUB,*999)
0730
0731 C Cluster jets with first hard parton
0732 DO I=1,NLJETS
0733 DO J=1,4
0734 PPM(J,I)=PJET(J,I)
0735 ENDDO
0736 ENDDO
0737
0738 NJETM=NLJETS
0739 DO IHARD=1,NLJETS
0740 NNM=NJETM+1
0741 DO J=1,4
0742 PPM(J,NNM)=p(J,IHARD)
0743 ENDDO
0744 CALL KTCLUS(KTSCHE,PPM,NNM,ECUT,YM,*999)
0745 IF(YM(NNM).GT.YCUT) THEN
0746 C Parton not clustered
0747 GOTO 90
0748 ENDIF
0749
0750 C Find jet clustered with parton
0751
0752 IP1=HIST(NNM)/NMAXKT
0753 IP2=MOD(HIST(NNM),NMAXKT)
0754 IF(IP2.NE.NNM.OR.IP1.LE.0)THEN
0755 GOTO 90
0756 ENDIF
0757 DO I=IP1,NJETM-1
0758 DO J=1,4
0759 PPM(J,I)=PPM(J,I+1)
0760 ENDDO
0761 ENDDO
0762 NJETM=NJETM-1
0763 ENDDO ! IHARD=1,NLJETS
0764 endif ! NLJETS.GT.MINJETS
0765
0766 DO I=1,MIN(NN,4)
0767 var2(1+I)=SQRT(Y(I))
0768 ENDDO
0769 C DKC WRITE(15,4001) (var2(I),I=1,nvar2)
0770
0771 90 CONTINUE
0772
0773 C Now perform jet clustering at the value chosen in qcut
0774
0775 CALL KTCLUS(KTSCHE,PP,NN,ECUT,Y,*999)
0776
0777 YCUT=qcut**2
0778 NCJET=0
0779
0780 C Reconstruct jet momenta
0781 CALL KTRECO(MOD(KTSCHE,10),PP,NN,ECUT,YCUT,YCUT,PJET,JET,
0782 $ NCJET,NSUB,*999)
0783
0784 ELSE IF (NN.EQ.1) THEN
0785
0786 Y(1)=PP(1,1)**2+PP(2,1)**2
0787 IF(Y(1).GT.YCUT)THEN
0788 NCJET=1
0789 DO I=1,4
0790 PJET(I,1)=PP(I,1)
0791 ENDDO
0792 ENDIF
0793 endif
0794
0795 if(idbg.eq.1) then
0796 write(LNHOUT,*) 'JETS'
0797 do i=1,ncjet
0798 PTJET =SQRT(PJET(1,i)**2+PJET(2,i)**2)
0799 ETAJET=PSERAP(PJET(1,i))
0800 PHIJET=ATAN2(PJET(2,i),PJET(1,i))
0801 write(LNHOUT,*) ptjet,etajet,phijet
0802 enddo
0803 endif
0804
0805 IF(NCJET.LT.NLJETS) THEN
0806 if(idbg.eq.1)
0807 $ WRITE(LNHOUT,*) 'Failed due to NCJET ',NCJET,' < ',NLJETS
0808 GOTO 999
0809 endif
0810
0811 C...Right number of jets - but the right jets?
0812 C For max. multiplicity case, count jets only to the NHARD:th jet
0813 IF(IEXC.EQ.0)THEN
0814 IF(NLJETS.GT.0)THEN
0815 YCUT=Y(NLJETS)
0816 CALL KTRECO(MOD(KTSCHE,10),PP,NN,ECUT,YCUT,YCUT,PJET,JET,
0817 $ NCJET,NSUB,*999)
0818 IF(clfact.GE.0d0) THEN
0819 CALL ALPSOR(PTCLUS,nljets,KPT,2)
0820 YCUT=MAX(qcut,PTCLUS(KPT(1)))**2
0821 ENDIF
0822 ENDIF
0823 ELSE IF(NCJET.GT.NLJETS) THEN
0824 if(idbg.eq.1)
0825 $ WRITE(LNHOUT,*) 'Failed due to NCJET ',NCJET,' > ',NLJETS
0826 GOTO 999
0827 ENDIF
0828
0829 C Cluster jets with hard partons, one at a time
0830 DO I=1,NLJETS
0831 DO J=1,4
0832 PPM(J,I)=PJET(J,I)
0833 ENDDO
0834 ENDDO
0835
0836 NJETM=NLJETS
0837 IF(clfact.NE.0) YCUT=clfact**2*YCUT
0838 c YCUT=qcut**2
0839 c YCUT=(1.5*qcut)**2
0840
0841 DO 120 IHARD=1,NLJETS
0842 NN=NJETM+1
0843 DO J=1,4
0844 PPM(J,NN)=p(J,IHARD)
0845 ENDDO
0846 CALL KTCLUS(KTSCHE,PPM,NN,ECUT,Y,*999)
0847
0848 IF(Y(NN).GT.YCUT) THEN
0849 C Parton not clustered
0850 if(idbg.eq.1)
0851 $ WRITE(LNHOUT,*) 'Failed due to parton ',IHARD,
0852 $ ' not clustered.'
0853 GOTO 999
0854 ENDIF
0855
0856 C Find jet clustered with parton
0857
0858 IP1=HIST(NN)/NMAXKT
0859 IP2=MOD(HIST(NN),NMAXKT)
0860 IF(IP2.NE.NN.OR.IP1.LE.0)THEN
0861 if(idbg.eq.1)
0862 $ WRITE(LNHOUT,*) 'Failed due to parton ',IHARD,
0863 $ ' not clustered.'
0864 GOTO 999
0865 ENDIF
0866 C Remove jet clustered with parton
0867 DO I=IP1,NJETM-1
0868 DO J=1,4
0869 PPM(J,I)=PPM(J,I+1)
0870 ENDDO
0871 ENDDO
0872 NJETM=NJETM-1
0873 120 CONTINUE
0874
0875 endif ! pt-ordered showers
0876 endif ! qcut.gt.0
0877
0878 C...Cluster particles with |eta| < etaclmax for histograms
0879 NN=0
0880 DO IHEP=1,NHEP
0881 IF (ISTHEP(IHEP).EQ.1
0882 & .AND.(ABS(IDHEP(IHEP)).LT.6.OR.IDHEP(IHEP).EQ.21)) THEN
0883 PTJET=sqrt(PHEP(1,IHEP)**2+PHEP(2,IHEP)**2)
0884 ETAJET=ABS(LOG(MIN((SQRT(PTJET**2+PHEP(3,IHEP)**2)+
0885 $ ABS(PHEP(3,IHEP)))/PTJET,1d5)))
0886 IF(ETAJET.GT.etaclmax) cycle
0887 NN=NN+1
0888 IF (NN.GT.NMAX) then
0889 CALL PYLIST(2)
0890 PRINT *, 'Too many particles: ', NN
0891 NN=NN-1
0892 GOTO 20
0893 ENDIF
0894 DO I=1,4
0895 PP(I,NN)=PHEP(I,IHEP)
0896 ENDDO
0897 ELSE if(idbg.eq.1)THEN
0898 PRINT *,'Skipping particle ',IHEP,ISTHEP(IHEP),IDHEP(IHEP)
0899 ENDIF
0900 ENDDO
0901
0902 20 ECUT=1
0903 IF (NN.GT.1) THEN
0904 CALL KTCLUS(KTSCHE,PP,NN,ECUT,Y,*999)
0905 ELSE IF(NN.EQ.1) THEN
0906 Y(1)=SQRT(PP(1,NN)**2+PP(2,NN)**2)
0907 ENDIF
0908
0909 DO I=1,MIN(NN,4)
0910 varev(46+I)=SQRT(Y(I))
0911 ENDDO
0912
0913
0914 RETURN
0915 4001 FORMAT(50E15.6)
0916 c HERWIG/PYTHIA TERMINATION:
0917 999 IPVETO=1
0918 END
0919
0920 C*********************************************************************
0921 C PYMASC
0922 C Implementation of scale used in Pythia parton showers
0923 C*********************************************************************
0924 SUBROUTINE PYMASC(scale)
0925 IMPLICIT NONE
0926
0927 C...Arguments
0928 REAL*8 scale
0929
0930 C...Functions
0931 REAL*8 SMDOT5
0932
0933 C...User process initialization commonblock.
0934 INTEGER MAXPUP
0935 PARAMETER (MAXPUP=100)
0936 INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
0937 DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
0938 COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
0939 & IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP),
0940 & LPRUP(MAXPUP)
0941 C...User process event common block.
0942 INTEGER MAXNUP
0943 PARAMETER (MAXNUP=500)
0944 INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
0945 DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
0946 COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP),
0947 & ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP),
0948 & VTIMUP(MAXNUP),SPINUP(MAXNUP)
0949
0950 C...Extra commonblock to transfer run info.
0951 INTEGER LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
0952 COMMON/UPPRIV/LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
0953
0954 C...Local variables
0955 INTEGER ICC1,ICC2,IJ,IDC1,IDC2,IC,IC1,IC2
0956 REAL*8 QMIN,QTMP
0957
0958 C Just use the scale read off the event record
0959 scale=SCALUP
0960
0961 C Alternatively:
0962
0963 C... Guesses for the correct scale
0964 C Assumptions:
0965 C (1) if the initial state is a color singlet, then
0966 C use s-hat for the scale
0967 C
0968 C (2) if color flow to the final state, use the minimum
0969 C of the dot products of color connected pairs
0970 C (times two for consistency with above)
0971
0972 QMIN=SMDOT5(PUP(1,1),PUP(1,2))
0973 ICC1=1
0974 ICC2=2
0975 C
0976 C For now, there is no generic way to guarantee the "right"
0977 C scale choice. Here, we take the HERWIG pt. of view and
0978 C choose the dot product of the colored connected "primary"
0979 C pairs.
0980 C
0981
0982 DO 101 IJ=1,NUP
0983 IF(MOTHUP(2,IJ).GT.2) GOTO 101
0984 IDC1=ICOLUP(1,IJ)
0985 IDC2=ICOLUP(2,IJ)
0986 IF(IDC1.EQ.0) IDC1=-1
0987 IF(IDC2.EQ.0) IDC2=-2
0988
0989 DO 201 IC=IJ+1,NUP
0990 IF(MOTHUP(2,IC).GT.2) GOTO 201
0991 IC1=ICOLUP(1,IC)
0992 IC2=ICOLUP(2,IC)
0993 IF(ISTUP(IC)*ISTUP(IJ).GE.1) THEN
0994 IF(IDC1.EQ.IC2.OR.IDC2.EQ.IC1) THEN
0995 QTMP=SMDOT5(PUP(1,IJ),PUP(1,IC))
0996 IF(QTMP.LT.QMIN) THEN
0997 QMIN=QTMP
0998 ICC1=IJ
0999 ICC2=IC
1000 ENDIF
1001 ENDIF
1002 ELSEIF(ISTUP(IC)*ISTUP(IJ).LE.-1) THEN
1003 IF(IDC1.EQ.IC1.OR.IDC2.EQ.IC2) THEN
1004 QTMP=SMDOT5(PUP(1,IJ),PUP(1,IC))
1005 IF(QTMP.LT.QMIN) THEN
1006 QMIN=QTMP
1007 ICC1=IJ
1008 ICC2=IC
1009 ENDIF
1010 ENDIF
1011 ENDIF
1012 201 CONTINUE
1013 101 CONTINUE
1014
1015 scale=QMIN
1016
1017 RETURN
1018 END
1019
1020 C...SMDOT5
1021 C Helper function
1022
1023 FUNCTION SMDOT5(V1,V2)
1024 IMPLICIT NONE
1025 REAL*8 SMDOT5,TEMP
1026 REAL*8 V1(5),V2(5)
1027 INTEGER I
1028
1029 SMDOT5=0D0
1030 TEMP=V1(4)*V2(4)
1031 DO I=1,3
1032 TEMP=TEMP-V1(I)*V2(I)
1033 ENDDO
1034
1035 SMDOT5=SQRT(ABS(TEMP))
1036
1037 RETURN
1038 END
1039
1040 C*********************************************************************
1041 C...set_matching
1042 C...Sets parameters for the matching, i.e. cuts and jet multiplicities
1043 C*********************************************************************
1044
1045 SUBROUTINE set_matching(npara,param,value)
1046 implicit none
1047 c
1048 c arguments
1049 c
1050 integer npara
1051 character*20 param(*),value(*)
1052
1053 C...User process initialization commonblock.
1054 INTEGER MAXPUP
1055 PARAMETER (MAXPUP=100)
1056 INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
1057 DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
1058 COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
1059 & IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP),
1060 & LPRUP(MAXPUP)
1061
1062 C...User process event common block.
1063 INTEGER MAXNUP
1064 PARAMETER (MAXNUP=500)
1065 INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
1066 DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
1067 COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP),
1068 & ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP),
1069 & VTIMUP(MAXNUP),SPINUP(MAXNUP)
1070
1071 C...Extra commonblock to transfer run info.
1072 INTEGER LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
1073 COMMON/UPPRIV/LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
1074
1075 C...Inputs for the matching algorithm
1076 double precision etcjet,rclmax,etaclmax,qcut,clfact
1077 integer maxjets,minjets,iexcfile,ktsche,nexcres,excres(30)
1078 common/MEMAIN/etcjet,rclmax,etaclmax,qcut,clfact,
1079 $ maxjets,minjets,iexcfile,ktsche,nexcres,excres
1080 CDKC DATA ktsche,maxjets,minjets,nexcres/0,-1,-1,0/
1081 DATA ktsche,maxjets,minjets,nexcres/0,-1,-1,0/
1082 DATA ktsche,nexcres/0,0/
1083 DATA qcut,clfact/0d0,0d0/
1084
1085 C...Commonblock to transfer event-by-event matching info
1086 INTEGER NLJETS,IEXC,Ifile
1087 DOUBLE PRECISION PTCLUS
1088 COMMON/MEMAEV/PTCLUS(20),NLJETS,IEXC,Ifile
1089
1090 C CALSIM VARIABLES
1091 INTEGER NCY,NCPHI
1092 DOUBLE PRECISION PI,DELY,DELPHI,ET,CTHCAL,STHCAL,CPHCAL,SPHCAL,
1093 $ YCMIN,YCMAX
1094 PARAMETER (NCY=50)
1095 PARAMETER (NCPHI=60,PI=3.141593D0)
1096 COMMON/CALOR/DELY,DELPHI,ET(NCY,NCPHI),
1097 $ CTHCAL(NCY),STHCAL(NCY),CPHCAL(NCPHI),SPHCAL(NCPHI),
1098 $ YCMIN,YCMAX
1099
1100 C...Local variables
1101 INTEGER I,MAXNJ,NREAD,MINJ,MAXJ
1102 parameter(MAXNJ=6)
1103 DOUBLE PRECISION XSTOT(MAXNJ),XSECTOT
1104 DOUBLE PRECISION ptjmin,etajmax,drjmin,ptbmin,etabmax,xqcut
1105
1106 CDKC -- initalize by assignment instead of data statement, whose behaviour is compiler dependent
1107 C minjets=-1
1108 C maxjets=-1
1109
1110 C...Need lower scale for final state radiation in e+e-
1111 IF(IABS(IDBMUP(1)).EQ.11.AND.IABS(IDBMUP(2)).EQ.11) then
1112 CALL PYGIVE('PARP(71)=1')
1113 ENDIF
1114
1115 C...CRUCIAL FOR JET-PARTON MATCHING: CALL UPVETO, ALLOW JET-PARTON MATCHING
1116 C call pygive('MSTP(143)=1')
1117
1118 if(ickkw.eq.1) then
1119
1120 C...Run PYPTFS instead of PYSHOW
1121 c CALL PYGIVE("MSTJ(41)=12")
1122
1123 c***********************************************************************
1124 c Read in jet cuts
1125 c***********************************************************************
1126
1127 call get_real (npara,param,value," ptj " ,ptjmin,7d3)
1128 call get_real (npara,param,value," etaj " ,etajmax,7d3)
1129 call get_real (npara,param,value," ptb " ,ptbmin,7d3)
1130 call get_real (npara,param,value," etab " ,etabmax,7d3)
1131 call get_real (npara,param,value," drjj " ,drjmin,7d3)
1132 call get_real (npara,param,value," xqcut " ,xqcut,0d0)
1133
1134 if(qcut.lt.xqcut) qcut=max(xqcut*1.2,xqcut+5)
1135 if(xqcut.le.0)then
1136 write(*,*) 'Warning! ME generation QCUT = 0. QCUT set to 0!'
1137 qcut=0
1138 endif
1139
1140 etajmax=min(etajmax,etabmax)
1141 ptjmin=max(ptjmin,ptbmin)
1142
1143 c IF(ICKKW.EQ.1) THEN
1144 c WRITE(*,*) ' '
1145 c WRITE(*,*) 'INPUT 0 FOR INCLUSIVE JET SAMPLE, 1 FOR EXCLUSIVE'
1146 c WRITE(*,*) '(SELECT 0 FOR HIGHEST PARTON MULTIPLICITY SAMPLE)'
1147 c WRITE(*,*) '(SELECT 1 OTHERWISE)'
1148 c READ(*,*) IEXCFILE
1149 c ENDIF
1150
1151 C INPUT PARAMETERS FOR CONE ALGORITHM
1152
1153 ETCJET=MAX(PTJMIN+5,1.2*PTJMIN)
1154 RCLMAX=DRJMIN
1155 ETACLMAX=ETAJMAX
1156 IF(qcut.le.0)THEN
1157 WRITE(*,*) 'JET CONE PARAMETERS FOR MATCHING:'
1158 WRITE(*,*) 'ET>',ETCJET,' R=',RCLMAX
1159 WRITE(*,*) 'DR(PARTON-JET)<',1.5*RCLMAX
1160 WRITE(*,*) 'ETA(JET)<',ETACLMAX
1161 ELSE
1162 WRITE(*,*) 'KT JET PARAMETERS FOR MATCHING:'
1163 WRITE(*,*) 'QCUT=',qcut
1164 WRITE(*,*) 'ETA(JET)<',ETACLMAX
1165 WRITE(*,*) 'Note that in ME generation, qcut = ',xqcut
1166 ENDIF
1167 endif
1168 return
1169 end
1170
1171 subroutine get_real(npara,param,value,name,var,def_value)
1172 c----------------------------------------------------------------------------------
1173 c finds the parameter named "name" in param and associate to "value" in value
1174 c----------------------------------------------------------------------------------
1175 implicit none
1176
1177 c
1178 c arguments
1179 c
1180 integer npara
1181 character*20 param(*),value(*)
1182 character*(*) name
1183 real*8 var,def_value
1184 c
1185 c local
1186 c
1187 logical found
1188 integer i
1189 c
1190 c start
1191 c
1192 i=1
1193 found=.false.
1194 do while(.not.found.and.i.le.npara)
1195 found = (index(param(i),name).ne.0)
1196 if (found) read(value(i),*) var
1197 c if (found) write (*,*) name,var
1198 i=i+1
1199 enddo
1200 if (.not.found) then
1201 write (*,*) "Warning: parameter ",name," not found"
1202 write (*,*) " setting it to default value ",def_value
1203 var=def_value
1204 else
1205 write(*,*),'Found parameter ',name,var
1206 endif
1207 return
1208
1209 end
1210 c
1211
1212 subroutine get_integer(npara,param,value,name,var,def_value)
1213 c----------------------------------------------------------------------------------
1214 c finds the parameter named "name" in param and associate to "value" in value
1215 c----------------------------------------------------------------------------------
1216 implicit none
1217 c
1218 c arguments
1219 c
1220 integer npara
1221 character*20 param(*),value(*)
1222 character*(*) name
1223 integer var,def_value
1224 c
1225 c local
1226 c
1227 logical found
1228 integer i
1229 c
1230 c start
1231 c
1232 i=1
1233 found=.false.
1234 do while(.not.found.and.i.le.npara)
1235 found = (index(param(i),name).ne.0)
1236 if (found) read(value(i),*) var
1237 c if (found) write (*,*) name,var
1238 i=i+1
1239 enddo
1240 if (.not.found) then
1241 write (*,*) "Warning: parameter ",name," not found"
1242 write (*,*) " setting it to default value ",def_value
1243 var=def_value
1244 else
1245 write(*,*)'Found parameter ',name,var
1246 endif
1247 return
1248
1249 end
1250
1251 C-----------------------------------------------------------------------
1252 SUBROUTINE ALPSOR(A,N,K,IOPT)
1253 C-----------------------------------------------------------------------
1254 C Sort A(N) into ascending order
1255 C IOPT = 1 : return sorted A and index array K
1256 C IOPT = 2 : return index array K only
1257 C-----------------------------------------------------------------------
1258 DOUBLE PRECISION A(N),B(5000)
1259 INTEGER N,I,J,IOPT,K(N),IL(5000),IR(5000)
1260 IF (N.GT.5000) then
1261 write(*,*) 'Too many entries to sort in alpsrt, stop'
1262 stop
1263 endif
1264 if(n.le.0) return
1265 IL(1)=0
1266 IR(1)=0
1267 DO 10 I=2,N
1268 IL(I)=0
1269 IR(I)=0
1270 J=1
1271 2 IF(A(I).GT.A(J)) GOTO 5
1272 3 IF(IL(J).EQ.0) GOTO 4
1273 J=IL(J)
1274 GOTO 2
1275 4 IR(I)=-J
1276 IL(J)=I
1277 GOTO 10
1278 5 IF(IR(J).LE.0) GOTO 6
1279 J=IR(J)
1280 GOTO 2
1281 6 IR(I)=IR(J)
1282 IR(J)=I
1283 10 CONTINUE
1284 I=1
1285 J=1
1286 GOTO 8
1287 20 J=IL(J)
1288 8 IF(IL(J).GT.0) GOTO 20
1289 9 K(I)=J
1290 B(I)=A(J)
1291 I=I+1
1292 IF(IR(J)) 12,30,13
1293 13 J=IR(J)
1294 GOTO 8
1295 12 J=-IR(J)
1296 GOTO 9
1297 30 IF(IOPT.EQ.2) RETURN
1298 DO 31 I=1,N
1299 31 A(I)=B(I)
1300 999 END
1301
1302 C-----------------------------------------------------------------------
1303 C----Calorimeter simulation obtained from Frank Paige 23 March 1988-----
1304 C
1305 C USE
1306 C
1307 C CALL CALINIM
1308 C CALL CALSIMM
1309 C
1310 C THEN TO FIND JETS WITH A SIMPLIFIED VERSION OF THE UA1 JET
1311 C ALGORITHM WITH JET RADIUS RJET AND MINIMUM SCALAR TRANSVERSE
1312 C ENERGY EJCUT
1313 C (RJET=1., EJCUT=5. FOR UA1)
1314 C USE
1315 C
1316 C CALL GETJETM(RJET,EJCUT)
1317 C
1318 C
1319 C-----------------------------------------------------------------------
1320 C
1321 C ADDED BY MIKE SEYMOUR: PARTON-LEVEL CALORIMETER. ALL PARTONS
1322 C ARE CONSIDERED TO BE HADRONS, SO IN FACT RESEM IS IGNORED
1323 C
1324 C CALL CALPARM
1325 C
1326 C HARD PARTICLE CALORIMETER. ONLY USES THOSE PARTICLES WHICH
1327 C CAME FROM THE HARD PROCESS, AND NOT THE UNDERLYING EVENT
1328 C
1329 C CALL CALHARM
1330 C
1331 C-----------------------------------------------------------------------
1332 SUBROUTINE CALINIM
1333 C
1334 C INITIALIZE CALORIMETER FOR CALSIMM AND GETJETM. NOTE THAT
1335 C BECAUSE THE INITIALIZATION IS SEPARATE, CALSIMM CAN BE
1336 C CALLED MORE THAN ONCE TO SIMULATE PILEUP OF SEVERAL EVENTS.
1337 C
1338 IMPLICIT NONE
1339 C...GETJET commonblocks
1340 INTEGER MNCY,MNCPHI,NCY,NCPHI,NJMAX,JETNO,NCJET
1341 DOUBLE PRECISION YCMIN,YCMAX,DELY,DELPHI,ET,STHCAL,CTHCAL,CPHCAL,
1342 & SPHCAL,PCJET,ETJET
1343 PARAMETER (MNCY=200)
1344 PARAMETER (MNCPHI=200)
1345 COMMON/CALORM/DELY,DELPHI,ET(MNCY,MNCPHI),
1346 $CTHCAL(MNCY),STHCAL(MNCY),CPHCAL(MNCPHI),SPHCAL(MNCPHI),
1347 $YCMIN,YCMAX,NCY,NCPHI
1348 PARAMETER (NJMAX=500)
1349 COMMON/GETCOMM/PCJET(4,NJMAX),ETJET(NJMAX),JETNO(MNCY,MNCPHI),
1350 $ NCJET
1351
1352 INTEGER IPHI,IY
1353 DOUBLE PRECISION PI,PHIX,YX,THX
1354 PARAMETER (PI=3.141593D0)
1355 LOGICAL FSTCAL
1356 DATA FSTCAL/.TRUE./
1357 C
1358 C INITIALIZE ET ARRAY.
1359 DO 100 IPHI=1,NCPHI
1360 DO 100 IY=1,NCY
1361 100 ET(IY,IPHI)=0.
1362 C
1363 IF (FSTCAL) THEN
1364 C CALCULATE TRIG. FUNCTIONS.
1365 DELPHI=2.*PI/FLOAT(NCPHI)
1366 DO 200 IPHI=1,NCPHI
1367 PHIX=DELPHI*(IPHI-.5)
1368 CPHCAL(IPHI)=COS(PHIX)
1369 SPHCAL(IPHI)=SIN(PHIX)
1370 200 CONTINUE
1371 DELY=(YCMAX-YCMIN)/FLOAT(NCY)
1372 DO 300 IY=1,NCY
1373 YX=DELY*(IY-.5)+YCMIN
1374 THX=2.*ATAN(EXP(-YX))
1375 CTHCAL(IY)=COS(THX)
1376 STHCAL(IY)=SIN(THX)
1377 300 CONTINUE
1378 FSTCAL=.FALSE.
1379 ENDIF
1380 END
1381 C
1382 SUBROUTINE CALSIMM
1383 C
1384 C SIMPLE CALORIMETER SIMULATION. ASSUME UNIFORM Y AND PHI
1385 C BINS
1386 C...HEPEVT commonblock.
1387 INTEGER NMXHEP,NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
1388 PARAMETER (NMXHEP=4000)
1389 COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
1390 &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
1391 DOUBLE PRECISION PHEP,VHEP
1392 SAVE /HEPEVT/
1393
1394 C...GETJET commonblocks
1395 INTEGER MNCY,MNCPHI,NCY,NCPHI,NJMAX,JETNO,NCJET
1396 DOUBLE PRECISION YCMIN,YCMAX,DELY,DELPHI,ET,STHCAL,CTHCAL,CPHCAL,
1397 & SPHCAL,PCJET,ETJET
1398 PARAMETER (MNCY=200)
1399 PARAMETER (MNCPHI=200)
1400 COMMON/CALORM/DELY,DELPHI,ET(MNCY,MNCPHI),
1401 $CTHCAL(MNCY),STHCAL(MNCY),CPHCAL(MNCPHI),SPHCAL(MNCPHI),
1402 $YCMIN,YCMAX,NCY,NCPHI
1403 PARAMETER (NJMAX=500)
1404 COMMON/GETCOMM/PCJET(4,NJMAX),ETJET(NJMAX),JETNO(MNCY,MNCPHI),
1405 $ NCJET
1406
1407 INTEGER IHEP,ID,IY,IPHI
1408 DOUBLE PRECISION PI,YIP,PSERAP,PHIIP,EIP
1409 PARAMETER (PI=3.141593D0)
1410 C
1411 C FILL CALORIMETER
1412 C
1413 DO 200 IHEP=1,NHEP
1414 IF (ISTHEP(IHEP).EQ.1) THEN
1415 YIP=PSERAP(PHEP(1,IHEP))
1416 IF(YIP.LT.YCMIN.OR.YIP.GT.YCMAX) GOTO 200
1417 ID=ABS(IDHEP(IHEP))
1418 C---EXCLUDE TOP QUARK, LEPTONS, PROMPT PHOTONS
1419 IF ((ID.GE.11.AND.ID.LE.16).OR.ID.EQ.6.OR.ID.EQ.22) GOTO 200
1420 C
1421 PHIIP=ATAN2(PHEP(2,IHEP),PHEP(1,IHEP))
1422 IF(PHIIP.LT.0.) PHIIP=PHIIP+2.*PI
1423 IY=INT((YIP-YCMIN)/DELY)+1
1424 IPHI=INT(PHIIP/DELPHI)+1
1425 EIP=PHEP(4,IHEP)
1426 C WEIGHT BY SIN(THETA)
1427 ET(IY,IPHI)=ET(IY,IPHI)+EIP*STHCAL(IY)
1428 ENDIF
1429 200 CONTINUE
1430 999 END
1431 SUBROUTINE GETJETM(RJET,EJCUT,ETAJCUT)
1432 C
1433 C SIMPLE JET-FINDING ALGORITHM (SIMILAR TO UA1).
1434 C
1435 C FIND HIGHEST REMAINING CELL > ETSTOP AND SUM SURROUNDING
1436 C CELLS WITH--
1437 C DELTA(Y)**2+DELTA(PHI)**2<RJET**2
1438 C ET>ECCUT.
1439 C KEEP JETS WITH ET>EJCUT AND ABS(ETA)<ETAJCUT
1440 C THE UA1 PARAMETERS ARE RJET=1.0 AND EJCUT=5.0
1441 C
1442 IMPLICIT NONE
1443 C...GETJET commonblocks
1444 INTEGER MNCY,MNCPHI,NCY,NCPHI,NJMAX,JETNO,NCJET
1445 DOUBLE PRECISION YCMIN,YCMAX,DELY,DELPHI,ET,STHCAL,CTHCAL,CPHCAL,
1446 & SPHCAL,PCJET,ETJET
1447 PARAMETER (MNCY=200)
1448 PARAMETER (MNCPHI=200)
1449 COMMON/CALORM/DELY,DELPHI,ET(MNCY,MNCPHI),
1450 $CTHCAL(MNCY),STHCAL(MNCY),CPHCAL(MNCPHI),SPHCAL(MNCPHI),
1451 $YCMIN,YCMAX,NCY,NCPHI
1452 PARAMETER (NJMAX=500)
1453 COMMON/GETCOMM/PCJET(4,NJMAX),ETJET(NJMAX),JETNO(MNCY,MNCPHI),
1454 $ NCJET
1455
1456 INTEGER IPHI,IY,J,K,NPHI1,NPHI2,NY1,
1457 & NY2,IPASS,IYMX,IPHIMX,ITLIS,IPHI1,IPHIX,IY1,IYX
1458 DOUBLE PRECISION PI,RJET,
1459 & ETMAX,ETSTOP,RR,ECCUT,PX,EJCUT
1460 PARAMETER (PI=3.141593D0)
1461 DOUBLE PRECISION ETAJCUT,PSERAP
1462 C
1463 C PARAMETERS
1464 DATA ECCUT/0.1D0/
1465 DATA ETSTOP/1.5D0/
1466 DATA ITLIS/6/
1467 C
1468 C INITIALIZE
1469 C
1470 DO 100 IPHI=1,NCPHI
1471 DO 100 IY=1,NCY
1472 100 JETNO(IY,IPHI)=0
1473 DO 110 J=1,NJMAX
1474 ETJET(J)=0.
1475 DO 110 K=1,4
1476 110 PCJET(K,J)=0.
1477 NCJET=0
1478 NPHI1=RJET/DELPHI
1479 NPHI2=2*NPHI1+1
1480 NY1=RJET/DELY
1481 NY2=2*NY1+1
1482 IPASS=0
1483 C
1484 C FIND HIGHEST CELL REMAINING
1485 C
1486 1 ETMAX=0.
1487 DO 200 IPHI=1,NCPHI
1488 DO 210 IY=1,NCY
1489 IF(ET(IY,IPHI).LT.ETMAX) GOTO 210
1490 IF(JETNO(IY,IPHI).NE.0) GOTO 210
1491 ETMAX=ET(IY,IPHI)
1492 IYMX=IY
1493 IPHIMX=IPHI
1494 210 CONTINUE
1495 200 CONTINUE
1496 IF(ETMAX.LT.ETSTOP) RETURN
1497 C
1498 C SUM CELLS
1499 C
1500 IPASS=IPASS+1
1501 IF(IPASS.GT.NCY*NCPHI) THEN
1502 WRITE(ITLIS,8888) IPASS
1503 8888 FORMAT(//' ERROR IN GETJETM...IPASS > ',I6)
1504 RETURN
1505 ENDIF
1506 NCJET=NCJET+1
1507 IF(NCJET.GT.NJMAX) THEN
1508 WRITE(ITLIS,9999) NCJET
1509 9999 FORMAT(//' ERROR IN GETJETM...NCJET > ',I5)
1510 RETURN
1511 ENDIF
1512 DO 300 IPHI1=1,NPHI2
1513 IPHIX=IPHIMX-NPHI1-1+IPHI1
1514 IF(IPHIX.LE.0) IPHIX=IPHIX+NCPHI
1515 IF(IPHIX.GT.NCPHI) IPHIX=IPHIX-NCPHI
1516 DO 310 IY1=1,NY2
1517 IYX=IYMX-NY1-1+IY1
1518 IF(IYX.LE.0) GOTO 310
1519 IF(IYX.GT.NCY) GOTO 310
1520 IF(JETNO(IYX,IPHIX).NE.0) GOTO 310
1521 RR=(DELY*(IY1-NY1-1))**2+(DELPHI*(IPHI1-NPHI1-1))**2
1522 IF(RR.GT.RJET**2) GOTO 310
1523 IF(ET(IYX,IPHIX).LT.ECCUT) GOTO 310
1524 PX=ET(IYX,IPHIX)/STHCAL(IYX)
1525 C ADD CELL TO JET
1526 PCJET(1,NCJET)=PCJET(1,NCJET)+PX*STHCAL(IYX)*CPHCAL(IPHIX)
1527 PCJET(2,NCJET)=PCJET(2,NCJET)+PX*STHCAL(IYX)*SPHCAL(IPHIX)
1528 PCJET(3,NCJET)=PCJET(3,NCJET)+PX*CTHCAL(IYX)
1529 PCJET(4,NCJET)=PCJET(4,NCJET)+PX
1530 ETJET(NCJET)=ETJET(NCJET)+ET(IYX,IPHIX)
1531 JETNO(IYX,IPHIX)=NCJET
1532 310 CONTINUE
1533 300 CONTINUE
1534 C
1535 C DISCARD JET IF ET < EJCUT.
1536 C
1537 IF(ETJET(NCJET).GT.EJCUT.AND.ABS(PSERAP(PCJET(1,NCJET))).LT
1538 $ .ETAJCUT) GOTO 1
1539 ETJET(NCJET)=0.
1540 DO 400 K=1,4
1541 400 PCJET(K,NCJET)=0.
1542 NCJET=NCJET-1
1543 GOTO 1
1544 END
1545 C-----------------------------------------------------------------------
1546 SUBROUTINE CALDELM(ISTLO,ISTHI)
1547 C LABEL ALL PARTICLES WITH STATUS BETWEEN ISTLO AND ISTHI (UNTIL A
1548 C PARTICLE WITH STATUS ISTOP IS FOUND) AS FINAL-STATE, CALL CALSIMM
1549 C AND THEN PUT LABELS BACK TO NORMAL
1550 C-----------------------------------------------------------------------
1551 IMPLICIT NONE
1552 INTEGER MAXNUP
1553 PARAMETER(MAXNUP=500)
1554 C...HEPEVT commonblock.
1555 INTEGER NMXHEP,NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
1556 PARAMETER (NMXHEP=4000)
1557 COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
1558 &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
1559 DOUBLE PRECISION PHEP,VHEP
1560 SAVE /HEPEVT/
1561 INTEGER ISTOLD(NMXHEP),IHEP,IST,ISTLO,ISTHI,ISTOP,IMO,icount
1562
1563
1564 CALL CALSIMM
1565 END
1566
