diff --git a/GeneratorInterface/LHEInterface/plugins/cms_hwwarn.F b/GeneratorInterface/LHEInterface/plugins/cms_hwwarn.F deleted file mode 100644 index 1b865f352040c..0000000000000 --- a/GeneratorInterface/LHEInterface/plugins/cms_hwwarn.F +++ /dev/null @@ -1,92 +0,0 @@ -CDECK ID>, HWWARN. -*CMZ :- -26/04/91 10.18.58 by Bryan Webber -*-- Author : Bryan Webber -C -C -C---------------------------------------------------------------------- -C modified by Fabian Stoeckli, March 2007 -C in order to catch run stop for CMSSW -C -C----------------------------------------------------------------------- - SUBROUTINE HWWARN(SUBRTN,ICODE) -C----------------------------------------------------------------------- -C DEALS WITH ERRORS DURING EXECUTION -C SUBRTN = NAME OF CALLING SUBROUTINE -C ICODE = ERROR CODE: - -1 NONFATAL, KILL EVENT & PRINT NOTHING -C 0- 49 NONFATAL, PRINT WARNING & CONTINUE -C 50- 99 NONFATAL, PRINT WARNING & JUMP -C 100-199 NONFATAL, DUMP & KILL EVENT -C 200-299 FATAL, TERMINATE RUN -C 300-399 FATAL, DUMP EVENT & TERMINATE RUN -C 400-499 FATAL, DUMP EVENT & STOP DEAD -C 500- FATAL, STOP DEAD WITH NO DUMP -C----------------------------------------------------------------------- - IMPLICIT NONE - DOUBLE PRECISION ZERO,ONE,TWO,THREE,FOUR,HALF - PARAMETER (ZERO =0.D0, ONE =1.D0, TWO =2.D0, - & THREE=3.D0, FOUR=4.D0, HALF=0.5D0) - -C - INTEGER NMXHEP,IERROR,NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP,IDHW, - & ISTAT,LWEVT,MAXER,MAXPR,NRN,NUMER,NUMERU,NWGTS - - DOUBLE PRECISION PHEP,VHEP,AVWGT,EVWGT,GAMWT,TLOUT,WBIGST,WGTMAX, - & WGTSUM,WSQSUM - - LOGICAL NOWGT,GENSOF - - PARAMETER (NMXHEP=4000) - COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP), - & JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP) - - COMMON/HWEVNT/AVWGT,EVWGT,GAMWT,TLOUT,WBIGST,WGTMAX,WGTSUM,WSQSUM, - & IDHW(NMXHEP),IERROR,ISTAT,LWEVT,MAXER,MAXPR,NOWGT,NRN(2),NUMER, - & NUMERU,NWGTS,GENSOF - - INTEGER ICODE - CHARACTER*6 SUBRTN - - IF (ICODE.GE.0) WRITE (6,10) SUBRTN,ICODE - 10 FORMAT(/' HWWARN CALLED FROM SUBPROGRAM ',A6,': CODE =',I4) - IF (ICODE.LT.0) THEN - IERROR=ICODE - RETURN - ELSEIF (ICODE.LT.100) THEN - WRITE (6,20) NEVHEP,NRN,EVWGT - 20 FORMAT(' EVENT',I8,': SEEDS =',I11,' &',I11, - &' WEIGHT =',E11.4/' EVENT SURVIVES. EXECUTION CONTINUES') - IF (ICODE.GT.49) RETURN - ELSEIF (ICODE.LT.200) THEN - WRITE (6,30) NEVHEP,NRN,EVWGT - 30 FORMAT(' EVENT',I8,': SEEDS =',I11,' &',I11, - &' WEIGHT =',E11.4/' EVENT KILLED. EXECUTION CONTINUES') - IERROR=ICODE - RETURN - ELSEIF (ICODE.LT.300) THEN - WRITE (6,40) - 40 FORMAT(' EVENT SURVIVES. RUN ENDS GRACEFULLY') - CALL HWEFIN - CALL HWAEND - CALL cmsending(ICODE) - ELSEIF (ICODE.LT.400) THEN - WRITE (6,50) - 50 FORMAT(' EVENT KILLED: DUMP FOLLOWS. RUN ENDS GRACEFULLY') - IERROR=ICODE - CALL HWUEPR - CALL HWUBPR - CALL HWEFIN - CALL HWAEND - CALL cmsending(ICODE) - ELSEIF (ICODE.LT.500) THEN - WRITE (6,60) - 60 FORMAT(' EVENT KILLED: DUMP FOLLOWS. RUN STOPS DEAD') - IERROR=ICODE - CALL HWUEPR - CALL HWUBPR - CALL cmsending(ICODE) - ELSE - WRITE (6,70) - 70 FORMAT(' RUN CANNOT CONTINUE') - CALL cmsending(ICODE) - ENDIF - END diff --git a/GeneratorInterface/LHEInterface/plugins/getjet.f b/GeneratorInterface/LHEInterface/plugins/getjet.f deleted file mode 100644 index aceddf32c2a3d..0000000000000 --- a/GeneratorInterface/LHEInterface/plugins/getjet.f +++ /dev/null @@ -1,279 +0,0 @@ -C----------------------------------------------------------------------- -C----Calorimeter simulation obtained from Frank Paige 23 March 1988----- -C -C USE -C -C CALL CALINI -C CALL CALSIM -C -C THEN TO FIND JETS WITH A SIMPLIFIED VERSION OF THE UA1 JET -C ALGORITHM WITH JET RADIUS RJET AND MINIMUM SCALAR TRANSVERSE -C ENERGY EJCUT -C (RJET=1., EJCUT=5. FOR UA1) -C USE -C -C CALL GETJET(RJET,EJCUT) -C -C -C----------------------------------------------------------------------- -C -C ADDED BY MIKE SEYMOUR: PARTON-LEVEL CALORIMETER. ALL PARTONS -C ARE CONSIDERED TO BE HADRONS, SO IN FACT RESEM IS IGNORED -C -C CALL CALPAR -C -C HARD PARTICLE CALORIMETER. ONLY USES THOSE PARTICLES WHICH -C CAME FROM THE HARD PROCESS, AND NOT THE UNDERLYING EVENT -C -C CALL CALHAR -C -C----------------------------------------------------------------------- - SUBROUTINE CALINI -C -C INITIALIZE CALORIMETER FOR CALSIM AND GETJET. NOTE THAT -C BECAUSE THE INITIALIZATION IS SEPARATE, CALSIM CAN BE -C CALLED MORE THAN ONCE TO SIMULATE PILEUP OF SEVERAL EVENTS. -C - IMPLICIT NONE -C...GETJET commonblocks - INTEGER MNCY,MNCPHI,NCY,NCPHI,NJMAX,JETNO,NCJET - DOUBLE PRECISION YCMIN,YCMAX,DELY,DELPHI,ET,STHCAL,CTHCAL,CPHCAL, - & SPHCAL,PCJET,ETJET - PARAMETER (MNCY=200) - PARAMETER (MNCPHI=200) - COMMON/CALOR/DELY,DELPHI,ET(MNCY,MNCPHI), - $CTHCAL(MNCY),STHCAL(MNCY),CPHCAL(MNCPHI),SPHCAL(MNCPHI), - $YCMIN,YCMAX,NCY,NCPHI - PARAMETER (NJMAX=500) - COMMON/GETCOM/PCJET(4,NJMAX),ETJET(NJMAX),JETNO(MNCY,MNCPHI),NCJET - - INTEGER IPHI,IY - DOUBLE PRECISION PI,PHIX,YX,THX - PARAMETER (PI=3.141593D0) - LOGICAL FSTCAL - DATA FSTCAL/.TRUE./ -C -C INITIALIZE ET ARRAY. - DO 100 IPHI=1,NCPHI - DO 100 IY=1,NCY -100 ET(IY,IPHI)=0. -C - IF (FSTCAL) THEN -C CALCULATE TRIG. FUNCTIONS. - DELPHI=2.*PI/FLOAT(NCPHI) - DO 200 IPHI=1,NCPHI - PHIX=DELPHI*(IPHI-.5) - CPHCAL(IPHI)=COS(PHIX) - SPHCAL(IPHI)=SIN(PHIX) -200 CONTINUE - DELY=(YCMAX-YCMIN)/FLOAT(NCY) - DO 300 IY=1,NCY - YX=DELY*(IY-.5)+YCMIN - THX=2.*ATAN(EXP(-YX)) - CTHCAL(IY)=COS(THX) - STHCAL(IY)=SIN(THX) -300 CONTINUE - FSTCAL=.FALSE. - ENDIF - END -C - SUBROUTINE CALSIM -C -C SIMPLE CALORIMETER SIMULATION. ASSUME UNIFORM Y AND PHI -C BINS -C...HEPEVT commonblock. - INTEGER NMXHEP,NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP - PARAMETER (NMXHEP=4000) - COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP), - &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP) - DOUBLE PRECISION PHEP,VHEP - SAVE /HEPEVT/ - -C...GETJET commonblocks - INTEGER MNCY,MNCPHI,NCY,NCPHI,NJMAX,JETNO,NCJET - DOUBLE PRECISION YCMIN,YCMAX,DELY,DELPHI,ET,STHCAL,CTHCAL,CPHCAL, - & SPHCAL,PCJET,ETJET - PARAMETER (MNCY=200) - PARAMETER (MNCPHI=200) - COMMON/CALOR/DELY,DELPHI,ET(MNCY,MNCPHI), - $CTHCAL(MNCY),STHCAL(MNCY),CPHCAL(MNCPHI),SPHCAL(MNCPHI), - $YCMIN,YCMAX,NCY,NCPHI - PARAMETER (NJMAX=500) - COMMON/GETCOM/PCJET(4,NJMAX),ETJET(NJMAX),JETNO(MNCY,MNCPHI),NCJET - - INTEGER IHEP,ID,IY,IPHI - DOUBLE PRECISION PI,YIP,PSERAP,PHIIP,EIP - PARAMETER (PI=3.141593D0) -C -C FILL CALORIMETER -C - DO 200 IHEP=1,NHEP - IF (ISTHEP(IHEP).EQ.1) THEN - YIP=PSERAP(PHEP(1,IHEP)) - IF(YIP.LT.YCMIN.OR.YIP.GT.YCMAX) GOTO 200 - ID=ABS(IDHEP(IHEP)) -C---EXCLUDE TOP QUARK, LEPTONS, PROMPT PHOTONS - IF ((ID.GE.11.AND.ID.LE.16).OR.ID.EQ.6.OR.ID.EQ.22) GOTO 200 -C - PHIIP=ATAN2(PHEP(2,IHEP),PHEP(1,IHEP)) - IF(PHIIP.LT.0.) PHIIP=PHIIP+2.*PI - IY=INT((YIP-YCMIN)/DELY)+1 - IPHI=INT(PHIIP/DELPHI)+1 - EIP=PHEP(4,IHEP) -C WEIGHT BY SIN(THETA) - ET(IY,IPHI)=ET(IY,IPHI)+EIP*STHCAL(IY) - ENDIF - 200 CONTINUE - 999 END - SUBROUTINE GETJET(RJET,EJCUT,ETAJCUT) -C -C SIMPLE JET-FINDING ALGORITHM (SIMILAR TO UA1). -C -C FIND HIGHEST REMAINING CELL > ETSTOP AND SUM SURROUNDING -C CELLS WITH-- -C DELTA(Y)**2+DELTA(PHI)**2ECCUT. -C KEEP JETS WITH ET>EJCUT AND ABS(ETA) ',I6) - RETURN - ENDIF - NCJET=NCJET+1 - IF(NCJET.GT.NJMAX) THEN - WRITE(ITLIS,9999) NCJET -9999 FORMAT(//' ERROR IN GETJET...NCJET > ',I5) - RETURN - ENDIF - DO 300 IPHI1=1,NPHI2 - IPHIX=IPHIMX-NPHI1-1+IPHI1 - IF(IPHIX.LE.0) IPHIX=IPHIX+NCPHI - IF(IPHIX.GT.NCPHI) IPHIX=IPHIX-NCPHI - DO 310 IY1=1,NY2 - IYX=IYMX-NY1-1+IY1 - IF(IYX.LE.0) GOTO 310 - IF(IYX.GT.NCY) GOTO 310 - IF(JETNO(IYX,IPHIX).NE.0) GOTO 310 - RR=(DELY*(IY1-NY1-1))**2+(DELPHI*(IPHI1-NPHI1-1))**2 - IF(RR.GT.RJET**2) GOTO 310 - IF(ET(IYX,IPHIX).LT.ECCUT) GOTO 310 - PX=ET(IYX,IPHIX)/STHCAL(IYX) -C ADD CELL TO JET - PCJET(1,NCJET)=PCJET(1,NCJET)+PX*STHCAL(IYX)*CPHCAL(IPHIX) - PCJET(2,NCJET)=PCJET(2,NCJET)+PX*STHCAL(IYX)*SPHCAL(IPHIX) - PCJET(3,NCJET)=PCJET(3,NCJET)+PX*CTHCAL(IYX) - PCJET(4,NCJET)=PCJET(4,NCJET)+PX - ETJET(NCJET)=ETJET(NCJET)+ET(IYX,IPHIX) - JETNO(IYX,IPHIX)=NCJET -310 CONTINUE -300 CONTINUE -C -C DISCARD JET IF ET < EJCUT. -C - IF(ETJET(NCJET).GT.EJCUT.AND.ABS(PSERAP(PCJET(1,NCJET))).LT - $ .ETAJCUT) GOTO 1 - ETJET(NCJET)=0. - DO 400 K=1,4 -400 PCJET(K,NCJET)=0. - NCJET=NCJET-1 - GOTO 1 - END -C----------------------------------------------------------------------- - SUBROUTINE CALDEL(ISTLO,ISTHI) -C LABEL ALL PARTICLES WITH STATUS BETWEEN ISTLO AND ISTHI (UNTIL A -C PARTICLE WITH STATUS ISTOP IS FOUND) AS FINAL-STATE, CALL CALSIM -C AND THEN PUT LABELS BACK TO NORMAL -C----------------------------------------------------------------------- - IMPLICIT NONE - INTEGER MAXNUP - PARAMETER(MAXNUP=500) -C...HEPEVT commonblock. - INTEGER NMXHEP,NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP - PARAMETER (NMXHEP=4000) - COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP), - &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP) - DOUBLE PRECISION PHEP,VHEP - SAVE /HEPEVT/ - INTEGER ISTOLD(NMXHEP),IHEP,IST,ISTLO,ISTHI,ISTOP,IMO,icount - -C...Inputs for the matching algorithm - double precision etcjet,rclmax,etaclmax,qcut,qfact - integer maxjets,minjets,iexcfile,ktsche - common/MEMAIN/etcjet,rclmax,etaclmax,qcut,qfact, - $ maxjets,minjets,iexcfile,ktsche - - - CALL CALSIM - END -C----------------------------------------------------------------------- - FUNCTION PSERAP(P) -C PSEUDO-RAPIDITY (-LOG TAN THETA/2) -C----------------------------------------------------------------------- - DOUBLE PRECISION PSERAP,P(3),PT,PL,TINY,THETA - PARAMETER (TINY=1D-3) - PT=SQRT(P(1)**2+P(2)**2)+TINY - PL=P(3) - THETA=ATAN2(PT,PL) - PSERAP=-LOG(TAN(0.5*THETA)) - END -C----------------------------------------------------------------------- diff --git a/GeneratorInterface/LHEInterface/plugins/hwmsct_dummy.F b/GeneratorInterface/LHEInterface/plugins/hwmsct_dummy.F deleted file mode 100644 index 86c6cd4ae27e2..0000000000000 --- a/GeneratorInterface/LHEInterface/plugins/hwmsct_dummy.F +++ /dev/null @@ -1,27 +0,0 @@ -CCC ADDED BY FABIAN STOECKLI -CCC CAUSE C++ DOESN'T LIKE LOGICAL !!!! - - - DOUBLE PRECISION FUNCTION HWMSCT_DUMMY(KERNTRR) - - DOUBLE PRECISION ABORE, KERNTRR - LOGICAL DOWORK - EXTERNAL HWMSCT - - CALL HWMSCT(DOWORK) - - ABORE=0.0d0 - - IF(DOWORK) THEN - ABORE = 1.0d0 - ELSE - ABORE = 0.0d0 - ENDIF - - -C write(*,*) ABORE, DOWORK - - HWMSCT_DUMMY = ABORE - - RETURN - END diff --git a/GeneratorInterface/LHEInterface/plugins/hwudat.F b/GeneratorInterface/LHEInterface/plugins/hwudat.F deleted file mode 100644 index 8c844887c73bd..0000000000000 --- a/GeneratorInterface/LHEInterface/plugins/hwudat.F +++ /dev/null @@ -1 +0,0 @@ -#include "GeneratorInterface/Herwig6Interface/src/hwudat.f" diff --git a/GeneratorInterface/LHEInterface/plugins/ktclusdble.f b/GeneratorInterface/LHEInterface/plugins/ktclusdble.f deleted file mode 100644 index 027f05acfdbf7..0000000000000 --- a/GeneratorInterface/LHEInterface/plugins/ktclusdble.f +++ /dev/null @@ -1,1567 +0,0 @@ -C----------------------------------------------------------------------- -C----------------------------------------------------------------------- -C----------------------------------------------------------------------- -C KTCLUS: written by Mike Seymour, July 1992. -C Last modified November 2000. -C Please send comments or suggestions to Mike.Seymour@rl.ac.uk -C -C This is a general-purpose kt clustering package. -C It can handle ee, ep and pp collisions. -C It is loosely based on the program of Siggi Bethke. -C -C The time taken (on a 10MIP machine) is (0.2microsec)*N**3 -C where N is the number of particles. -C Over 90 percent of this time is used in subroutine KTPMIN, which -C simply finds the minimum member of a one-dimensional array. -C It is well worth thinking about optimization: on the SPARCstation -C a factor of two increase was obtained simply by increasing the -C optimization level from its default value. -C -C The approach is to separate the different stages of analysis. -C KTCLUS does all the clustering and records a merging history. -C It returns a simple list of the y values at which each merging -C occured. Then the following routines can be called to give extra -C information on the most recently analysed event. -C KTCLUR is identical but includes an R parameter, see below. -C KTYCUT gives the number of jets at each given YCUT value. -C KTYSUB gives the number of sub-jets at each given YCUT value. -C KTBEAM gives same info as KTCLUS but only for merges with the beam -C KTJOIN gives same info as KTCLUS but for merges of sub-jets. -C KTRECO reconstructs the jet momenta at a given value of YCUT. -C It also gives information on which jets at scale YCUT belong to -C which macro-jets at scale YMAC, for studying sub-jet properties. -C KTINCL reconstructs the jet momenta according to the inclusive jet -C definition of Ellis and Soper. -C KTISUB, KTIJOI and KTIREC are like KTYSUB, KTJOIN and KTRECO, -C except that they only apply to one inclusive jet at a time, -C with the pt of that jet automatically used for ECUT. -C KTWICH gives a list of which particles ended up in which jets. -C KTWCHS gives the same thing, but only for subjets. -C Note that the numbering of jets used by these two routines is -C guaranteed to be the same as that used by KTRECO. -C -C The collision type and analysis type are indicated by the first -C argument of KTCLUS. IMODE= where -C TYPE: 1=>ee, 2=>ep with p in -z direction, 3=>pe, 4=>pp -C ANGLE: 1=>angular kt def., 2=>DeltaR, 3=>f(DeltaEta,DeltaPhi) -C where f()=2(cosh(eta)-cos(phi)) is the QCD emission metric -C MONO: 1=>derive relative pseudoparticle angles from jets -C 2=>monotonic definitions of relative angles -C RECOM: 1=>E recombination scheme, 2=>pt scheme, 3=>pt**2 scheme -C -C There are also abbreviated forms for the most common combinations: -C IMODE=1 => E scheme in e+e- (=1111) -C 2 => E scheme in ep (=2111) -C 3 => E scheme in pe (=3111) -C 4 => E scheme in pp (=4111) -C 5 => covariant E scheme in pp (=4211) -C 6 => covariant pt-scheme in pp (=4212) -C 7 => covariant monotonic pt**2-scheme in pp (=4223) -C -C KTRECO no longer needs to reconstruct the momenta according to the -C same recombination scheme in which they were clustered. Its first -C argument gives the scheme, taking the same values as RECOM above. -C -C Note that unlike previous versions, all variables which hold y -C values have been named in a consistent way: -C Y() is the output scale at which jets were merged, -C YCUT is the input scale at which jets should be counted, and -C jet-momenta reconstructed etc, -C YMAC is the input macro-jet scale, used in determining whether -C or not each jet is a sub-jet. -C The original scheme defined in our papers is equivalent to always -C setting YMAC=1. -C Whenever a YCUT or YMAC variable is used, it is rounded down -C infinitesimally, so that for example, setting YCUT=Y(2) refers -C to the scale where the event is 2-jet, even if rounding errors -C have shifted its value slightly. -C -C An R parameter can be used in hadron-hadron collisions by -C calling KTCLUR instead of KTCLUS. This is as suggested by -C Ellis and Soper, but implemented slightly differently, -C as in M.H. Seymour, LU TP 94/2 (submitted to Nucl. Phys. B.). -C R**2 multiplies the single Kt everywhere it is used. -C Calling KTCLUR with R=1 is identical to calling KTCLUS. -C R plays a similar role to the jet radius in a cone-type algorithm, -C but is scaled up by about 40% (ie R=0.7 in a cone algorithm is -C similar to this algorithm with R=1). -C Note that R.EQ.1 must be used for the e+e- and ep versions, -C and is strongly recommended for the hadron-hadron version. -C However, R values smaller than 1 have been found to be useful for -C certain applications, particularly the mass reconstruction of -C highly-boosted colour-singlets such as high-pt hadronic Ws, -C as in M.H. Seymour, LU TP 93/8 (to appear in Z. Phys. C.). -C Situations in which R<1 is useful are likely to also be those in -C which the inclusive reconstruction method is more useful. -C -C Also included is a set of routines for doing Lorentz boosts: -C KTLBST finds the boost matrix to/from the cm frame of a 4-vector -C KTRROT finds the rotation matrix from one vector to another -C KTMMUL multiplies together two matrices -C KTVMUL multiplies a vector by a matrix -C KTINVT inverts a transformation matrix (nb NOT a general 4 by 4) -C KTFRAM boosts a list of vectors between two arbitrary frames -C KTBREI boosts a list of vectors between the lab and Breit frames -C KTHADR boosts a list of vectors between the lab and hadronic cmf -C The last two need the momenta in the +z direction of the lepton -C and hadron beams, and the 4-momentum of the outgoing lepton. -C -C The main reference is: -C S. Catani, Yu.L. Dokshitzer, M.H. Seymour and B.R. Webber, -C Nucl.Phys.B406(1993)187. -C The ep version was proposed in: -C S. Catani, Yu.L. Dokshitzer and B.R. Webber, -C Phys.Lett.285B(1992)291. -C The inclusive reconstruction method was proposed in: -C S.D. Ellis and D.E. Soper, -C Phys.Rev.D48(1993)3160. -C -C----------------------------------------------------------------------- -C----------------------------------------------------------------------- -C----------------------------------------------------------------------- - SUBROUTINE KTCLUS(IMODE,PP,NN,ECUT,Y,*) - IMPLICIT NONE -C---DO CLUSTER ANALYSIS OF PARTICLES IN PP -C -C IMODE = INPUT : DESCRIBED ABOVE -C PP(I,J) = INPUT : 4-MOMENTUM OF Jth PARTICLE: I=1,4 => PX,PY,PZ,E -C NN = INPUT : NUMBER OF PARTICLES -C ECUT = INPUT : DENOMINATOR OF KT MEASURE. IF ZERO, ETOT IS USED -C Y(J) = OUTPUT : VALUE OF Y FOR WHICH EVENT CHANGES FROM BEING -C J JET TO J-1 JET -C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT -C COULD NOT BE PROCESSED (MOST LIKELY DUE TO TOO MANY PARTICLES) -C -C NOTE THAT THE MOMENTA ARE DECLARED DOUBLE PRECISION, -C AND ALL OTHER FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION -C - INTEGER IMODE,NN - DOUBLE PRECISION PP(4,*) - DOUBLE PRECISION ECUT,Y(*),ONE - ONE=1 - CALL KTCLUR(IMODE,PP,NN,ONE,ECUT,Y,*999) - RETURN - 999 RETURN 1 - END -C----------------------------------------------------------------------- - SUBROUTINE KTCLUR(IMODE,PP,NN,R,ECUT,Y,*) - IMPLICIT NONE -C---DO CLUSTER ANALYSIS OF PARTICLES IN PP -C -C IMODE = INPUT : DESCRIBED ABOVE -C PP(I,J) = INPUT : 4-MOMENTUM OF Jth PARTICLE: I=1,4 => PX,PY,PZ,E -C NN = INPUT : NUMBER OF PARTICLES -C R = INPUT : ELLIS AND SOPER'S R PARAMETER, SEE ABOVE. -C ECUT = INPUT : DENOMINATOR OF KT MEASURE. IF ZERO, ETOT IS USED -C Y(J) = OUTPUT : VALUE OF Y FOR WHICH EVENT CHANGES FROM BEING -C J JET TO J-1 JET -C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT -C COULD NOT BE PROCESSED (MOST LIKELY DUE TO TOO MANY PARTICLES) -C -C NOTE THAT THE MOMENTA ARE DECLARED DOUBLE PRECISION, -C AND ALL OTHER FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION -C - INTEGER NMAX,IM,IMODE,TYPE,ANGL,MONO,RECO,N,I,J,NN, - & IMIN,JMIN,KMIN,NUM,HIST,INJET,IABBR,NABBR - PARAMETER (NMAX=512,NABBR=7) - DOUBLE PRECISION PP(4,*) - DOUBLE PRECISION R,ECUT,Y(*),P,KT,ETOT,RSQ,KTP,KTS,KTPAIR,KTSING, - & KTMIN,ETSQ,KTLAST,KTMAX,KTTMP - LOGICAL FIRST - CHARACTER TITLE(4,4)*10 -C---KT RECORDS THE KT**2 OF EACH MERGING. -C---KTLAST RECORDS FOR EACH MERGING, THE HIGHEST ECUT**2 FOR WHICH THE -C RESULT IS NOT MERGED WITH THE BEAM (COULD BE LARGER THAN THE -C KT**2 AT WHICH IT WAS MERGED IF THE KT VALUES ARE NOT MONOTONIC). -C THIS MAY SOUND POINTLESS, BUT ITS USEFUL FOR DETERMINING WHETHER -C SUB-JETS SURVIVED TO SCALE Y=YMAC OR NOT. -C---HIST RECORDS MERGING HISTORY: -C N=>DELETED TRACK N, M*NMAX+N=>MERGED TRACKS M AND N (M PX,PY,PZ,E -C NN = INPUT : NUMBER OF PARTICLES -C ECUT = INPUT : DENOMINATOR OF KT MEASURE. IF ZERO, ETOT IS USED -C YCUT = INPUT : Y VALUE AT WHICH TO RECONSTRUCT JET MOMENTA -C YMAC = INPUT : Y VALUE USED TO DEFINE MACRO-JETS, TO DETERMINE -C WHICH JETS ARE SUB-JETS -C PJET(I,J)=OUTPUT : 4-MOMENTUM OF Jth JET AT SCALE YCUT -C JET(J) =OUTPUT : THE MACRO-JET WHICH CONTAINS THE Jth JET, -C SET TO ZERO IF JET IS NOT A SUB-JET -C NJET =OUTPUT : THE NUMBER OF JETS -C NSUB =OUTPUT : THE NUMBER OF SUB-JETS (EQUAL TO THE NUMBER OF -C NON-ZERO ENTRIES IN JET()) -C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT -C COULD NOT BE PROCESSED -C -C NOTE THAT THE MOMENTA ARE DECLARED DOUBLE PRECISION, -C AND ALL OTHER FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION -C - INTEGER NMAX,RECO,NUM,N,NN,NJET,NSUB,JET(*),HIST,IMIN,JMIN,I,J - PARAMETER (NMAX=512) - DOUBLE PRECISION PP(4,*),PJET(4,*) - DOUBLE PRECISION ECUT,P,KT,KTP,KTS,ETOT,RSQ,ETSQ,YCUT,YMAC,KTLAST, - & ROUND - PARAMETER (ROUND=0.99999D0) - COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX), - & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM -C---CHECK INPUT - IF (RECO.LT.1.OR.RECO.GT.3) THEN - PRINT *,'RECO=',RECO - CALL KTWARN('KTRECO',100,*999) - ENDIF -C---COPY PP TO P - N=NN - IF (NUM.NE.NN) CALL KTWARN('KTRECO',101,*999) - CALL KTCOPY(PP,N,P,(RECO.NE.1)) - IF (ECUT.EQ.0) THEN - ETSQ=1/ETOT**2 - ELSE - ETSQ=1/ECUT**2 - ENDIF -C---KEEP MERGING UNTIL YCUT - 100 IF (ETSQ*KT(N).LT.ROUND*YCUT) THEN - IF (HIST(N).LE.NMAX) THEN - CALL KTMOVE(P,KTP,KTS,NMAX,N,HIST(N),0) - ELSE - IMIN=HIST(N)/NMAX - JMIN=HIST(N)-IMIN*NMAX - CALL KTMERG(P,KTP,KTS,NMAX,IMIN,JMIN,N,0,0,0,RECO) - CALL KTMOVE(P,KTP,KTS,NMAX,N,JMIN,0) - ENDIF - N=N-1 - IF (N.GT.0) GOTO 100 - ENDIF -C---IF YCUT IS TOO LARGE THERE ARE NO JETS - NJET=N - NSUB=N - IF (N.EQ.0) RETURN -C---SET UP OUTPUT MOMENTA - DO 210 I=1,NJET - IF (RECO.EQ.1) THEN - DO 200 J=1,4 - PJET(J,I)=P(J,I) - 200 CONTINUE - ELSE - PJET(1,I)=P(6,I)*COS(P(8,I)) - PJET(2,I)=P(6,I)*SIN(P(8,I)) - PJET(3,I)=P(6,I)*SINH(P(7,I)) - PJET(4,I)=P(6,I)*COSH(P(7,I)) - ENDIF - JET(I)=I - 210 CONTINUE -C---KEEP MERGING UNTIL YMAC TO FIND THE FATE OF EACH JET - 300 IF (ETSQ*KT(N).LT.ROUND*YMAC) THEN - IF (HIST(N).LE.NMAX) THEN - IMIN=0 - JMIN=HIST(N) - NSUB=NSUB-1 - ELSE - IMIN=HIST(N)/NMAX - JMIN=HIST(N)-IMIN*NMAX - IF (ETSQ*KTLAST(N).LT.ROUND*YMAC) NSUB=NSUB-1 - ENDIF - DO 310 I=1,NJET - IF (JET(I).EQ.JMIN) JET(I)=IMIN - IF (JET(I).EQ.N) JET(I)=JMIN - 310 CONTINUE - N=N-1 - IF (N.GT.0) GOTO 300 - ENDIF - RETURN - 999 RETURN 1 - END -C----------------------------------------------------------------------- - SUBROUTINE KTINCL(RECO,PP,NN,PJET,JET,NJET,*) - IMPLICIT NONE -C---RECONSTRUCT KINEMATICS OF JET SYSTEM, WHICH HAS ALREADY BEEN -C ANALYSED BY KTCLUS ACCORDING TO THE INCLUSIVE JET DEFINITION. NOTE -C THAT NO CONSISTENCY CHECK IS MADE: USER IS TRUSTED TO USE THE SAME -C PP VALUES AS FOR KTCLUS -C -C RECO = INPUT : RECOMBINATION SCHEME (NEED NOT BE SAME AS KTCLUS) -C PP(I,J) = INPUT : 4-MOMENTUM OF Jth PARTICLE: I=1,4 => PX,PY,PZ,E -C NN = INPUT : NUMBER OF PARTICLES -C PJET(I,J)=OUTPUT : 4-MOMENTUM OF Jth JET AT SCALE YCUT -C JET(J) =OUTPUT : THE JET WHICH CONTAINS THE Jth PARTICLE -C NJET =OUTPUT : THE NUMBER OF JETS -C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT -C COULD NOT BE PROCESSED -C -C NOTE THAT THE MOMENTA ARE DECLARED DOUBLE PRECISION, -C AND ALL OTHER FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION -C - INTEGER NMAX,RECO,NUM,N,NN,NJET,JET(*),HIST,IMIN,JMIN,I,J - PARAMETER (NMAX=512) - DOUBLE PRECISION PP(4,*),PJET(4,*) - DOUBLE PRECISION P,KT,KTP,KTS,ETOT,RSQ,KTLAST - COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX), - & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM -C---CHECK INPUT - IF (RECO.LT.1.OR.RECO.GT.3) CALL KTWARN('KTINCL',100,*999) -C---COPY PP TO P - N=NN - IF (NUM.NE.NN) CALL KTWARN('KTINCL',101,*999) - CALL KTCOPY(PP,N,P,(RECO.NE.1)) -C---INITIALLY EVERY PARTICLE IS IN ITS OWN JET - DO 100 I=1,NN - JET(I)=I - 100 CONTINUE -C---KEEP MERGING TO THE BITTER END - NJET=0 - 200 IF (N.GT.0) THEN - IF (HIST(N).LE.NMAX) THEN - IMIN=0 - JMIN=HIST(N) - NJET=NJET+1 - IF (RECO.EQ.1) THEN - DO 300 J=1,4 - PJET(J,NJET)=P(J,JMIN) - 300 CONTINUE - ELSE - PJET(1,NJET)=P(6,JMIN)*COS(P(8,JMIN)) - PJET(2,NJET)=P(6,JMIN)*SIN(P(8,JMIN)) - PJET(3,NJET)=P(6,JMIN)*SINH(P(7,JMIN)) - PJET(4,NJET)=P(6,JMIN)*COSH(P(7,JMIN)) - ENDIF - CALL KTMOVE(P,KTP,KTS,NMAX,N,JMIN,0) - ELSE - IMIN=HIST(N)/NMAX - JMIN=HIST(N)-IMIN*NMAX - CALL KTMERG(P,KTP,KTS,NMAX,IMIN,JMIN,N,0,0,0,RECO) - CALL KTMOVE(P,KTP,KTS,NMAX,N,JMIN,0) - ENDIF - DO 400 I=1,NN - IF (JET(I).EQ.JMIN) JET(I)=IMIN - IF (JET(I).EQ.N) JET(I)=JMIN - IF (JET(I).EQ.0) JET(I)=-NJET - 400 CONTINUE - N=N-1 - GOTO 200 - ENDIF -C---FINALLY EVERY PARTICLE MUST BE IN AN INCLUSIVE JET - DO 500 I=1,NN -C---IF THERE ARE ANY UNASSIGNED PARTICLES SOMETHING MUST HAVE GONE WRONG - IF (JET(I).GE.0) CALL KTWARN('KTINCL',102,*999) - JET(I)=-JET(I) - 500 CONTINUE - RETURN - 999 RETURN 1 - END -C----------------------------------------------------------------------- - SUBROUTINE KTISUB(N,NY,YCUT,NSUB,*) - IMPLICIT NONE -C---COUNT THE NUMBER OF SUB-JETS IN THE Nth INCLUSIVE JET OF AN EVENT -C THAT HAS ALREADY BEEN ANALYSED BY KTCLUS. -C -C N = INPUT : WHICH INCLUSIVE JET TO USE -C NY = INPUT : NUMBER OF YCUT VALUES -C YCUT(J) = INPUT : Y VALUES AT WHICH NUMBERS OF SUB-JETS ARE COUNTED -C NSUB(J) =OUTPUT : NUMBER OF SUB-JETS AT YCUT(J) -C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT -C COULD NOT BE PROCESSED -C -C NOTE THAT ALL FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION -C - INTEGER N,NY,NSUB(NY),NMAX,HIST,I,J,NUM,NM - PARAMETER (NMAX=512) - DOUBLE PRECISION YCUT(NY),ETOT,RSQ,P,KT,KTP,KTS,KTLAST,ROUND,EPS - PARAMETER (ROUND=0.99999D0) - COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX), - & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM - DATA EPS/1D-6/ - DO 100 I=1,NY - NSUB(I)=0 - 100 CONTINUE -C---FIND WHICH MERGING CORRESPONDS TO THE NTH INCLUSIVE JET - NM=0 - J=0 - DO 110 I=NUM,1,-1 - IF (HIST(I).LE.NMAX) J=J+1 - IF (J.EQ.N) THEN - NM=I - GOTO 120 - ENDIF - 110 CONTINUE - 120 CONTINUE -C---GIVE UP IF THERE ARE LESS THAN N INCLUSIVE JETS - IF (NM.EQ.0) CALL KTWARN('KTISUB',100,*999) - DO 210 I=NUM,1,-1 - DO 200 J=1,NY - IF (NSUB(J).EQ.0.AND.RSQ*KT(I).GE.ROUND*YCUT(J)*KT(NM)) - & NSUB(J)=I - IF (NSUB(J).NE.0.AND.ABS(KTLAST(I)-KTLAST(NM)).GT.EPS) - & NSUB(J)=NSUB(J)-1 - 200 CONTINUE - 210 CONTINUE - RETURN - 999 RETURN 1 - END -C----------------------------------------------------------------------- - SUBROUTINE KTIJOI(N,Y,*) - IMPLICIT NONE -C---GIVE SAME INFORMATION AS LAST CALL TO KTCLUS EXCEPT THAT ONLY -C MERGES OF TWO SUB-JETS INSIDE THE Nth INCLUSIVE JET ARE RECORDED -C -C N = INPUT : WHICH INCLUSIVE JET TO USE -C Y(J) =OUTPUT : Y VALUE WHERE JET CHANGED FROM HAVING -C J+1 SUB-JETS TO HAVING J -C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT -C COULD NOT BE PROCESSED -C -C NOTE THAT ALL FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION -C - INTEGER NMAX,HIST,NUM,I,J,N,NM - PARAMETER (NMAX=512) - DOUBLE PRECISION ETOT,RSQ,P,KT,KTP,KTS,Y(*),KTLAST,EPS - COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX), - & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM - DATA EPS/1D-6/ -C---FIND WHICH MERGING CORRESPONDS TO THE NTH INCLUSIVE JET - NM=0 - J=0 - DO 100 I=NUM,1,-1 - IF (HIST(I).LE.NMAX) J=J+1 - IF (J.EQ.N) THEN - NM=I - GOTO 105 - ENDIF - 100 CONTINUE - 105 CONTINUE -C---GIVE UP IF THERE ARE LESS THAN N INCLUSIVE JETS - IF (NM.EQ.0) CALL KTWARN('KTIJOI',100,*999) - J=1 - DO 110 I=1,NUM - IF (HIST(I).GT.NMAX.AND.ABS(KTLAST(I)-KTLAST(NM)).LT.EPS) THEN - Y(J)=RSQ*KT(I)/KT(NM) - J=J+1 - ENDIF - 110 CONTINUE - DO 200 I=J,NUM - Y(I)=0 - 200 CONTINUE - RETURN - 999 RETURN 1 - END -C----------------------------------------------------------------------- - SUBROUTINE KTIREC(RECO,PP,NN,N,YCUT,PSUB,NSUB,*) - IMPLICIT NONE -C---RECONSTRUCT KINEMATICS OF SUB-JET SYSTEM IN THE Nth INCLUSIVE JET -C OF AN EVENT THAT HAS ALREADY BEEN ANALYSED BY KTCLUS -C -C RECO = INPUT : RECOMBINATION SCHEME (NEED NOT BE SAME AS KTCLUS) -C PP(I,J) = INPUT : 4-MOMENTUM OF Jth PARTICLE: I=1,4 => PX,PY,PZ,E -C NN = INPUT : NUMBER OF PARTICLES -C N = INPUT : WHICH INCLUSIVE JET TO USE -C YCUT = INPUT : Y VALUE AT WHICH TO RECONSTRUCT JET MOMENTA -C PSUB(I,J)=OUTPUT : 4-MOMENTUM OF Jth SUB-JET AT SCALE YCUT -C NSUB =OUTPUT : THE NUMBER OF SUB-JETS -C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT -C COULD NOT BE PROCESSED -C -C NOTE THAT THE MOMENTA ARE DECLARED DOUBLE PRECISION, -C AND ALL OTHER FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION -C - INTEGER NMAX,RECO,NUM,NN,NJET,NSUB,JET,HIST,I,J,N,NM - PARAMETER (NMAX=512) - DOUBLE PRECISION PP(4,*),PSUB(4,*) - DOUBLE PRECISION ECUT,P,KT,KTP,KTS,ETOT,RSQ,YCUT,YMAC,KTLAST - COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX), - & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM - DIMENSION JET(NMAX) -C---FIND WHICH MERGING CORRESPONDS TO THE NTH INCLUSIVE JET - NM=0 - J=0 - DO 100 I=NUM,1,-1 - IF (HIST(I).LE.NMAX) J=J+1 - IF (J.EQ.N) THEN - NM=I - GOTO 110 - ENDIF - 100 CONTINUE - 110 CONTINUE -C---GIVE UP IF THERE ARE LESS THAN N INCLUSIVE JETS - IF (NM.EQ.0) CALL KTWARN('KTIREC',102,*999) -C---RECONSTRUCT THE JETS AT THE APPROPRIATE SCALE - ECUT=SQRT(KT(NM)/RSQ) - YMAC=RSQ - CALL KTRECO(RECO,PP,NN,ECUT,YCUT,YMAC,PSUB,JET,NJET,NSUB,*999) -C---GET RID OF THE ONES THAT DO NOT END UP IN THE JET WE WANT - NSUB=0 - DO 210 I=1,NJET - IF (JET(I).EQ.HIST(NM)) THEN - NSUB=NSUB+1 - DO 200 J=1,4 - PSUB(J,NSUB)=PSUB(J,I) - 200 CONTINUE - ENDIF - 210 CONTINUE - RETURN - 999 RETURN 1 - END -C----------------------------------------------------------------------- - SUBROUTINE KTWICH(ECUT,YCUT,JET,NJET,*) - IMPLICIT NONE -C---GIVE A LIST OF WHICH JET EACH ORIGINAL PARTICLE ENDED UP IN AT SCALE -C YCUT, TOGETHER WITH THE NUMBER OF JETS AT THAT SCALE. -C -C ECUT = INPUT : DENOMINATOR OF KT MEASURE. IF ZERO, ETOT IS USED -C YCUT = INPUT : Y VALUE AT WHICH TO DEFINE JETS -C JET(J) =OUTPUT : THE JET WHICH CONTAINS THE Jth PARTICLE, -C SET TO ZERO IF IT WAS PUT INTO THE BEAM JETS -C NJET =OUTPUT : THE NUMBER OF JETS AT SCALE YCUT (SO JET() -C ENTRIES WILL BE IN THE RANGE 0 -> NJET) -C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT -C COULD NOT BE PROCESSED -C -C NOTE THAT ALL FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION -C - INTEGER JET(*),NJET,NTEMP - DOUBLE PRECISION ECUT,YCUT - CALL KTWCHS(ECUT,YCUT,YCUT,JET,NJET,NTEMP,*999) - RETURN - 999 RETURN 1 - END -C----------------------------------------------------------------------- - SUBROUTINE KTWCHS(ECUT,YCUT,YMAC,JET,NJET,NSUB,*) - IMPLICIT NONE -C---GIVE A LIST OF WHICH SUB-JET EACH ORIGINAL PARTICLE ENDED UP IN AT -C SCALE YCUT, WITH MACRO-JET SCALE YMAC, TOGETHER WITH THE NUMBER OF -C JETS AT SCALE YCUT AND THE NUMBER OF THEM WHICH ARE SUB-JETS. -C -C ECUT = INPUT : DENOMINATOR OF KT MEASURE. IF ZERO, ETOT IS USED -C YCUT = INPUT : Y VALUE AT WHICH TO DEFINE JETS -C YMAC = INPUT : Y VALUE AT WHICH TO DEFINE MACRO-JETS -C JET(J) =OUTPUT : THE JET WHICH CONTAINS THE Jth PARTICLE, -C SET TO ZERO IF IT WAS PUT INTO THE BEAM JETS -C NJET =OUTPUT : THE NUMBER OF JETS AT SCALE YCUT (SO JET() -C ENTRIES WILL BE IN THE RANGE 0 -> NJET) -C NSUB =OUTPUT : THE NUMBER OF SUB-JETS AT SCALE YCUT, WITH -C MACRO-JETS DEFINED AT SCALE YMAC (SO ONLY NSUB -C OF THE JETS 1 -> NJET WILL APPEAR IN JET()) -C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT -C COULD NOT BE PROCESSED -C -C NOTE THAT ALL FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION -C - INTEGER NMAX,JET(*),NJET,NSUB,HIST,NUM,I,J,JSUB - PARAMETER (NMAX=512) - DOUBLE PRECISION P1(4,NMAX),P2(4,NMAX) - DOUBLE PRECISION ECUT,YCUT,YMAC,ZERO,ETOT,RSQ,P,KTP,KTS,KT,KTLAST - COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX), - & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM - DIMENSION JSUB(NMAX) -C---THE MOMENTA HAVE TO BEEN GIVEN LEGAL VALUES, -C EVEN THOUGH THEY WILL NEVER BE USED - DATA ((P1(J,I),I=1,NMAX),J=1,4),ZERO - & /NMAX*1,NMAX*0,NMAX*0,NMAX*1,0/ -C---FIRST GET A LIST OF WHICH PARTICLE IS IN WHICH JET AT YCUT - CALL KTRECO(1,P1,NUM,ECUT,ZERO,YCUT,P2,JET,NJET,NSUB,*999) -C---THEN FIND OUT WHICH JETS ARE SUBJETS - CALL KTRECO(1,P1,NUM,ECUT,YCUT,YMAC,P2,JSUB,NJET,NSUB,*999) -C---AND MODIFY JET() ACCORDINGLY - DO 10 I=1,NUM - IF (JET(I).NE.0) THEN - IF (JSUB(JET(I)).EQ.0) JET(I)=0 - ENDIF - 10 CONTINUE - RETURN - 999 RETURN 1 - END -C----------------------------------------------------------------------- - SUBROUTINE KTFRAM(IOPT,CMF,SIGN,Z,XZ,N,P,Q,*) - IMPLICIT NONE -C---BOOST PARTICLES IN P TO/FROM FRAME GIVEN BY CMF, Z, XZ. -C---IN THIS FRAME CMZ IS STATIONARY, -C Z IS ALONG THE (SIGN)Z-AXIS (SIGN=+ OR -) -C XZ IS IN THE X-Z PLANE (WITH POSITIVE X COMPONENT) -C---IF Z HAS LENGTH ZERO, OR SIGN=0, NO ROTATION IS PERFORMED -C---IF XZ HAS ZERO COMPONENT PERPENDICULAR TO Z IN THAT FRAME, -C NO AZIMUTHAL ROTATION IS PERFORMED -C -C IOPT = INPUT : 0=TO FRAME, 1=FROM FRAME -C CMF(I) = INPUT : 4-MOMENTUM WHICH IS STATIONARY IN THE FRAME -C SIGN = INPUT : DIRECTION OF Z IN THE FRAME, NOTE THAT -C ONLY ITS SIGN IS USED, NOT ITS MAGNITUDE -C Z(I) = INPUT : 4-MOMENTUM WHICH LIES ON THE (SIGN)Z-AXIS -C XZ(I) = INPUT : 4-MOMENTUM WHICH LIES IN THE X-Z PLANE -C N = INPUT : NUMBER OF PARTICLES IN P -C P(I,J) = INPUT : 4-MOMENTUM OF JTH PARTICLE BEFORE -C Q(I,J) = OUTPUT : 4-MOMENTUM OF JTH PARTICLE AFTER -C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT -C COULD NOT BE PROCESSED -C -C NOTE THAT ALL MOMENTA ARE DOUBLE PRECISION -C -C NOTE THAT IT IS SAFE TO CALL WITH P=Q -C - INTEGER IOPT,I,N - DOUBLE PRECISION CMF(4),SIGN,Z(4),XZ(4),P(4,N),Q(4,N), - & R(4,4),NEW(4),OLD(4) - IF (IOPT.LT.0.OR.IOPT.GT.1) CALL KTWARN('KTFRAM',200,*999) -C---FIND BOOST TO GET THERE FROM LAB - CALL KTUNIT(R) - CALL KTLBST(0,R,CMF,*999) -C---FIND ROTATION TO PUT BOOSTED Z ON THE (SIGN)Z AXIS - IF (SIGN.NE.0) THEN - CALL KTVMUL(R,Z,OLD) - IF (OLD(1).NE.0.OR.OLD(2).NE.0.OR.OLD(3).NE.0) THEN - NEW(1)=0 - NEW(2)=0 - NEW(3)=SIGN - NEW(4)=ABS(SIGN) - CALL KTRROT(R,OLD,NEW,*999) -C---FIND ROTATION TO PUT BOOSTED AND ROTATED XZ INTO X-Z PLANE - CALL KTVMUL(R,XZ,OLD) - IF (OLD(1).NE.0.OR.OLD(2).NE.0) THEN - NEW(1)=1 - NEW(2)=0 - NEW(3)=0 - NEW(4)=1 - OLD(3)=0 -C---NOTE THAT A POTENTIALLY AWKWARD SPECIAL CASE IS AVERTED, BECAUSE IF -C OLD AND NEW ARE EXACTLY BACK-TO-BACK, THE ROTATION AXIS IS UNDEFINED -C BUT IN THAT CASE KTRROT WILL USE THE Z AXIS, AS REQUIRED - CALL KTRROT(R,OLD,NEW,*999) - ENDIF - ENDIF - ENDIF -C---INVERT THE TRANSFORMATION IF NECESSARY - IF (IOPT.EQ.1) CALL KTINVT(R,R) -C---APPLY THE RESULT TO ALL THE VECTORS - DO 30 I=1,N - CALL KTVMUL(R,P(1,I),Q(1,I)) - 30 CONTINUE - RETURN - 999 RETURN 1 - END -C----------------------------------------------------------------------- - SUBROUTINE KTBREI(IOPT,PLEP,PHAD,POUT,N,P,Q,*) - IMPLICIT NONE -C---BOOST PARTICLES IN P TO/FROM BREIT FRAME -C -C IOPT = INPUT : 0/2=TO BREIT FRAME, 1/3=FROM BREIT FRAME -C 0/1=NO AZIMUTHAL ROTATION AFTERWARDS -C 2/3=LEPTON PLANE ROTATED INTO THE X-Z PLANE -C PLEP = INPUT : MOMENTUM OF INCOMING LEPTON IN +Z DIRECTION -C PHAD = INPUT : MOMENTUM OF INCOMING HADRON IN +Z DIRECTION -C POUT(I) = INPUT : 4-MOMENTUM OF OUTGOING LEPTON -C N = INPUT : NUMBER OF PARTICLES IN P -C P(I,J) = INPUT : 4-MOMENTUM OF JTH PARTICLE BEFORE -C Q(I,J) = OUTPUT : 4-MOMENTUM OF JTH PARTICLE AFTER -C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT -C COULD NOT BE PROCESSED (MOST LIKELY DUE TO PARTICLES HAVING SMALLER -C ENERGY THAN MOMENTUM) -C -C NOTE THAT ALL MOMENTA ARE DOUBLE PRECISION -C -C NOTE THAT IT IS SAFE TO CALL WITH P=Q -C - INTEGER IOPT,N - DOUBLE PRECISION PLEP,PHAD,POUT(4),P(4,N),Q(4,N), - & CMF(4),Z(4),XZ(4),DOT,QDQ -C---CHECK INPUT - IF (IOPT.LT.0.OR.IOPT.GT.3) CALL KTWARN('KTBREI',200,*999) -C---FIND 4-MOMENTUM OF BREIT FRAME (TIMES AN ARBITRARY FACTOR) - DOT=ABS(PHAD)*(ABS(PLEP)-POUT(4))-PHAD*(PLEP-POUT(3)) - QDQ=(ABS(PLEP)-POUT(4))**2-(PLEP-POUT(3))**2-POUT(2)**2-POUT(1)**2 - CMF(1)=DOT*( -POUT(1)) - CMF(2)=DOT*( -POUT(2)) - CMF(3)=DOT*( PLEP -POUT(3))-QDQ* PHAD - CMF(4)=DOT*(ABS(PLEP)-POUT(4))-QDQ*ABS(PHAD) -C---FIND ROTATION TO PUT INCOMING HADRON BACK ON Z-AXIS - Z(1)=0 - Z(2)=0 - Z(3)=PHAD - Z(4)=ABS(PHAD) - XZ(1)=0 - XZ(2)=0 - XZ(3)=0 - XZ(4)=0 -C---DO THE BOOST - IF (IOPT.LE.1) THEN - CALL KTFRAM(IOPT,CMF,PHAD,Z,XZ,N,P,Q,*999) - ELSE - CALL KTFRAM(IOPT-2,CMF,PHAD,Z,POUT,N,P,Q,*999) - ENDIF - RETURN - 999 RETURN 1 - END -C----------------------------------------------------------------------- - SUBROUTINE KTHADR(IOPT,PLEP,PHAD,POUT,N,P,Q,*) - IMPLICIT NONE -C---BOOST PARTICLES IN P TO/FROM HADRONIC CMF -C -C ARGUMENTS ARE EXACTLY AS FOR KTBREI -C -C NOTE THAT ALL MOMENTA ARE DOUBLE PRECISION -C -C NOTE THAT IT IS SAFE TO CALL WITH P=Q -C - INTEGER IOPT,N - DOUBLE PRECISION PLEP,PHAD,POUT(4),P(4,N),Q(4,N), - & CMF(4),Z(4),XZ(4) -C---CHECK INPUT - IF (IOPT.LT.0.OR.IOPT.GT.3) CALL KTWARN('KTHADR',200,*999) -C---FIND 4-MOMENTUM OF HADRONIC CMF - CMF(1)= -POUT(1) - CMF(2)= -POUT(2) - CMF(3)= PLEP -POUT(3)+ PHAD - CMF(4)=ABS(PLEP)-POUT(4)+ABS(PHAD) -C---FIND ROTATION TO PUT INCOMING HADRON BACK ON Z-AXIS - Z(1)=0 - Z(2)=0 - Z(3)=PHAD - Z(4)=ABS(PHAD) - XZ(1)=0 - XZ(2)=0 - XZ(3)=0 - XZ(4)=0 -C---DO THE BOOST - IF (IOPT.LE.1) THEN - CALL KTFRAM(IOPT,CMF,PHAD,Z,XZ,N,P,Q,*999) - ELSE - CALL KTFRAM(IOPT-2,CMF,PHAD,Z,POUT,N,P,Q,*999) - ENDIF - RETURN - 999 RETURN 1 - END -C----------------------------------------------------------------------- - FUNCTION KTPAIR(ANGL,P,Q,ANGLE) - IMPLICIT NONE -C---CALCULATE LOCAL KT OF PAIR, USING ANGULAR SCHEME: -C 1=>ANGULAR, 2=>DeltaR, 3=>f(DeltaEta,DeltaPhi) -C WHERE f(eta,phi)=2(COSH(eta)-COS(phi)) IS THE QCD EMISSION METRIC -C---IF ANGLE<0, IT IS SET TO THE ANGULAR PART OF THE LOCAL KT ON RETURN -C IF ANGLE>0, IT IS USED INSTEAD OF THE ANGULAR PART OF THE LOCAL KT - INTEGER ANGL,I - DOUBLE PRECISION P(9),Q(9),KTPAIR,R,KTMDPI,ANGLE,ETA,PHI,ESQ -C---COMPONENTS OF MOMENTA ARE PX,PY,PZ,E,1/P,PT,ETA,PHI,PT**2 - R=ANGLE - IF (ANGL.EQ.1) THEN - IF (R.LE.0) R=2*(1-(P(1)*Q(1)+P(2)*Q(2)+P(3)*Q(3))*(P(5)*Q(5))) - ESQ=MIN(P(4),Q(4))**2 - ELSEIF (ANGL.EQ.2.OR.ANGL.EQ.3) THEN - IF (R.LE.0) THEN - ETA=P(7)-Q(7) - PHI=KTMDPI(P(8)-Q(8)) - IF (ANGL.EQ.2) THEN - R=ETA**2+PHI**2 - ELSE - R=2*(COSH(ETA)-COS(PHI)) - ENDIF - ENDIF - ESQ=MIN(P(9),Q(9)) - ELSEIF (ANGL.EQ.4) THEN - ESQ=(1d0/(P(5)*Q(5))-P(1)*Q(1)-P(2)*Q(2)- - &P(3)*Q(3))*2D0/(P(5)*Q(5))/(0.0001D0+1d0/P(5)+1d0/Q(5))**2 - R=1d0 - ELSE - CALL KTWARN('KTPAIR',200,*999) - STOP - ENDIF - KTPAIR=ESQ*R - IF (ANGLE.LT.0) ANGLE=R - 999 END -C----------------------------------------------------------------------- - FUNCTION KTSING(ANGL,TYPE,P) - IMPLICIT NONE -C---CALCULATE KT OF PARTICLE, USING ANGULAR SCHEME: -C 1=>ANGULAR, 2=>DeltaR, 3=>f(DeltaEta,DeltaPhi) -C---TYPE=1 FOR E+E-, 2 FOR EP, 3 FOR PE, 4 FOR PP -C FOR EP, PROTON DIRECTION IS DEFINED AS -Z -C FOR PE, PROTON DIRECTION IS DEFINED AS +Z - INTEGER ANGL,TYPE - DOUBLE PRECISION P(9),KTSING,COSTH,R,SMALL - DATA SMALL/1D-4/ - IF (ANGL.EQ.1.OR.ANGL.EQ.4) THEN - COSTH=P(3)*P(5) - IF (TYPE.EQ.2) THEN - COSTH=-COSTH - ELSEIF (TYPE.EQ.4) THEN - COSTH=ABS(COSTH) - ELSEIF (TYPE.NE.1.AND.TYPE.NE.3) THEN - CALL KTWARN('KTSING',200,*999) - STOP - ENDIF - R=2*(1-COSTH) -C---IF CLOSE TO BEAM, USE APPROX 2*(1-COS(THETA))=SIN**2(THETA) - IF (R.LT.SMALL) R=(P(1)**2+P(2)**2)*P(5)**2 - KTSING=P(4)**2*R - ELSEIF (ANGL.EQ.2.OR.ANGL.EQ.3) THEN - KTSING=P(9) - ELSE - CALL KTWARN('KTSING',201,*999) - STOP - ENDIF - 999 END -C----------------------------------------------------------------------- - SUBROUTINE KTPMIN(A,NMAX,N,IMIN,JMIN) - IMPLICIT NONE -C---FIND THE MINIMUM MEMBER OF A(NMAX,NMAX) WITH IMIN < JMIN <= N - INTEGER NMAX,N,IMIN,JMIN,KMIN,I,J,K -C---REMEMBER THAT A(X+(Y-1)*NMAX)=A(X,Y) -C THESE LOOPING VARIABLES ARE J=Y-2, I=X+(Y-1)*NMAX - DOUBLE PRECISION A(*),AMIN - K=1+NMAX - KMIN=K - AMIN=A(KMIN) - DO 110 J=0,N-2 - DO 100 I=K,K+J - IF (A(I).LT.AMIN) THEN - KMIN=I - AMIN=A(KMIN) - ENDIF - 100 CONTINUE - K=K+NMAX - 110 CONTINUE - JMIN=KMIN/NMAX+1 - IMIN=KMIN-(JMIN-1)*NMAX - END -C----------------------------------------------------------------------- - SUBROUTINE KTSMIN(A,NMAX,N,IMIN) - IMPLICIT NONE -C---FIND THE MINIMUM MEMBER OF A - INTEGER N,NMAX,IMIN,I - DOUBLE PRECISION A(NMAX) - IMIN=1 - DO 100 I=1,N - IF (A(I).LT.A(IMIN)) IMIN=I - 100 CONTINUE - END -C----------------------------------------------------------------------- - SUBROUTINE KTCOPY(A,N,B,ONSHLL) - IMPLICIT NONE -C---COPY FROM A TO B. 5TH=1/(3-MTM), 6TH=PT, 7TH=ETA, 8TH=PHI, 9TH=PT**2 -C IF ONSHLL IS .TRUE. PARTICLE ENTRIES ARE PUT ON-SHELL BY SETTING E=P - INTEGER I,N - DOUBLE PRECISION A(4,N) - LOGICAL ONSHLL - DOUBLE PRECISION B(9,N),ETAMAX,SINMIN,EPS - DATA ETAMAX,SINMIN,EPS/10,0,1D-6/ -C---SINMIN GETS CALCULATED ON FIRST CALL - IF (SINMIN.EQ.0) SINMIN=1/COSH(ETAMAX) - DO 100 I=1,N - B(1,I)=A(1,I) - B(2,I)=A(2,I) - B(3,I)=A(3,I) - B(4,I)=A(4,I) - B(5,I)=SQRT(A(1,I)**2+A(2,I)**2+A(3,I)**2) - IF (ONSHLL) B(4,I)=B(5,I) - IF (B(5,I).EQ.0) B(5,I)=1D-10 - B(5,I)=1/B(5,I) - B(9,I)=A(1,I)**2+A(2,I)**2 - B(6,I)=SQRT(B(9,I)) - B(7,I)=B(6,I)*B(5,I) - IF (B(7,I).GT.SINMIN) THEN - B(7,I)=A(4,I)**2-A(3,I)**2 - IF (B(7,I).LE.EPS*B(4,I)**2.OR.ONSHLL) B(7,I)=B(9,I) - B(7,I)=LOG((B(4,I)+ABS(B(3,I)))**2/B(7,I))/2 - ELSE - B(7,I)=ETAMAX+2 - ENDIF - B(7,I)=SIGN(B(7,I),B(3,I)) - IF (A(1,I).EQ.0 .AND. A(2,I).EQ.0) THEN - B(8,I)=0 - ELSE - B(8,I)=ATAN2(A(2,I),A(1,I)) - ENDIF - 100 CONTINUE - END -C----------------------------------------------------------------------- - SUBROUTINE KTMERG(P,KTP,KTS,NMAX,I,J,N,TYPE,ANGL,MONO,RECO) - IMPLICIT NONE -C---MERGE THE Jth PARTICLE IN P INTO THE Ith PARTICLE -C J IS ASSUMED GREATER THAN I. P CONTAINS N PARTICLES BEFORE MERGING. -C---ALSO RECALCULATING THE CORRESPONDING KTP AND KTS VALUES IF MONO.GT.0 -C FROM THE RECOMBINED ANGULAR MEASURES IF MONO.GT.1 -C---NOTE THAT IF MONO.LE.0, TYPE AND ANGL ARE NOT USED - INTEGER ANGL,RECO,TYPE,I,J,K,N,NMAX,MONO - DOUBLE PRECISION P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),PT,PTT, - & KTMDPI,KTUP,PI,PJ,ANG,KTPAIR,KTSING,ETAMAX,EPS - KTUP(I,J)=KTP(MAX(I,J),MIN(I,J)) - DATA ETAMAX,EPS/10,1D-6/ - IF (J.LE.I) CALL KTWARN('KTMERG',200,*999) -C---COMBINE ANGULAR MEASURES IF NECESSARY - IF (MONO.GT.1) THEN - DO 100 K=1,N - IF (K.NE.I.AND.K.NE.J) THEN - IF (RECO.EQ.1) THEN - PI=P(4,I) - PJ=P(4,J) - ELSEIF (RECO.EQ.2) THEN - PI=P(6,I) - PJ=P(6,J) - ELSEIF (RECO.EQ.3) THEN - PI=P(9,I) - PJ=P(9,J) - ELSE - CALL KTWARN('KTMERG',201,*999) - STOP - ENDIF - IF (PI.EQ.0.AND.PJ.EQ.0) THEN - PI=1 - PJ=1 - ENDIF - KTP(MAX(I,K),MIN(I,K))= - & (PI*KTUP(I,K)+PJ*KTUP(J,K))/(PI+PJ) - ENDIF - 100 CONTINUE - ENDIF - IF (RECO.EQ.1) THEN -C---VECTOR ADDITION - P(1,I)=P(1,I)+P(1,J) - P(2,I)=P(2,I)+P(2,J) - P(3,I)=P(3,I)+P(3,J) -c P(4,I)=P(4,I)+P(4,J) ! JA - P(5,I)=SQRT(P(1,I)**2+P(2,I)**2+P(3,I)**2) - P(4,I)=P(5,I) ! JA (Massless scheme) - IF (P(5,I).EQ.0) THEN - P(5,I)=1 - ELSE - P(5,I)=1/P(5,I) - ENDIF - ELSEIF (RECO.EQ.2) THEN -C---PT WEIGHTED ETA-PHI ADDITION - PT=P(6,I)+P(6,J) - IF (PT.EQ.0) THEN - PTT=1 - ELSE - PTT=1/PT - ENDIF - P(7,I)=(P(6,I)*P(7,I)+P(6,J)*P(7,J))*PTT - P(8,I)=KTMDPI(P(8,I)+P(6,J)*PTT*KTMDPI(P(8,J)-P(8,I))) - P(6,I)=PT - P(9,I)=PT**2 - ELSEIF (RECO.EQ.3) THEN -C---PT**2 WEIGHTED ETA-PHI ADDITION - PT=P(9,I)+P(9,J) - IF (PT.EQ.0) THEN - PTT=1 - ELSE - PTT=1/PT - ENDIF - P(7,I)=(P(9,I)*P(7,I)+P(9,J)*P(7,J))*PTT - P(8,I)=KTMDPI(P(8,I)+P(9,J)*PTT*KTMDPI(P(8,J)-P(8,I))) - P(6,I)=P(6,I)+P(6,J) - P(9,I)=P(6,I)**2 - ELSE - CALL KTWARN('KTMERG',202,*999) - STOP - ENDIF -C---IF MONO.GT.0 CALCULATE NEW KT MEASURES. IF MONO.GT.1 USE ANGULAR ONES. - IF (MONO.LE.0) RETURN -C---CONVERTING BETWEEN 4-MTM AND PT,ETA,PHI IF NECESSARY - IF (ANGL.NE.1.AND.RECO.EQ.1) THEN - P(9,I)=P(1,I)**2+P(2,I)**2 - P(7,I)=P(4,I)**2-P(3,I)**2 - IF (P(7,I).LE.EPS*P(4,I)**2) P(7,I)=P(9,I) - IF (P(7,I).GT.0) THEN - P(7,I)=LOG((P(4,I)+ABS(P(3,I)))**2/P(7,I))/2 - IF (P(7,I).GT.ETAMAX) P(7,I)=ETAMAX+2 - ELSE - P(7,I)=ETAMAX+2 - ENDIF - P(7,I)=SIGN(P(7,I),P(3,I)) - IF (P(1,I).NE.0.AND.P(2,I).NE.0) THEN - P(8,I)=ATAN2(P(2,I),P(1,I)) - ELSE - P(8,I)=0 - ENDIF - ELSEIF (ANGL.EQ.1.AND.RECO.NE.1) THEN - P(1,I)=P(6,I)*COS(P(8,I)) - P(2,I)=P(6,I)*SIN(P(8,I)) - P(3,I)=P(6,I)*SINH(P(7,I)) - P(4,I)=P(6,I)*COSH(P(7,I)) - IF (P(4,I).NE.0) THEN - P(5,I)=1/P(4,I) - ELSE - P(5,I)=1 - ENDIF - ENDIF - ANG=0 - DO 200 K=1,N - IF (K.NE.I.AND.K.NE.J) THEN - IF (MONO.GT.1) ANG=KTUP(I,K) - KTP(MIN(I,K),MAX(I,K))= - & KTPAIR(ANGL,P(1,I),P(1,K),ANG) - ENDIF - 200 CONTINUE - KTS(I)=KTSING(ANGL,TYPE,P(1,I)) - 999 END -C----------------------------------------------------------------------- - SUBROUTINE KTMOVE(P,KTP,KTS,NMAX,N,J,IOPT) - IMPLICIT NONE -C---MOVE THE Nth PARTICLE IN P TO THE Jth POSITION -C---ALSO MOVING KTP AND KTS IF IOPT.GT.0 - INTEGER I,J,N,NMAX,IOPT - DOUBLE PRECISION P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX) - DO 100 I=1,9 - P(I,J)=P(I,N) - 100 CONTINUE - IF (IOPT.LE.0) RETURN - DO 110 I=1,J-1 - KTP(I,J)=KTP(I,N) - KTP(J,I)=KTP(N,I) - 110 CONTINUE - DO 120 I=J+1,N-1 - KTP(J,I)=KTP(I,N) - KTP(I,J)=KTP(N,I) - 120 CONTINUE - KTS(J)=KTS(N) - END -C----------------------------------------------------------------------- - SUBROUTINE KTUNIT(R) - IMPLICIT NONE -C SET R EQUAL TO THE 4 BY 4 IDENTITY MATRIX - DOUBLE PRECISION R(4,4) - INTEGER I,J - DO 20 I=1,4 - DO 10 J=1,4 - R(I,J)=0 - IF (I.EQ.J) R(I,J)=1 - 10 CONTINUE - 20 CONTINUE - END -C----------------------------------------------------------------------- - SUBROUTINE KTLBST(IOPT,R,A,*) - IMPLICIT NONE -C PREMULTIPLY R BY THE 4 BY 4 MATRIX TO -C LORENTZ BOOST TO/FROM THE CM FRAME OF A -C IOPT=0 => TO -C IOPT=1 => FROM -C -C LAST ARGUMENT IS LABEL TO JUMP TO IF A IS NOT TIME-LIKE -C - INTEGER IOPT,I,J - DOUBLE PRECISION R(4,4),A(4),B(4),C(4,4),M - DO 10 I=1,4 - B(I)=A(I) - 10 CONTINUE - M=B(4)**2-B(1)**2-B(2)**2-B(3)**2 - IF (M.LE.0) CALL KTWARN('KTLBST',100,*999) - M=SQRT(M) - B(4)=B(4)+M - M=1/(M*B(4)) - IF (IOPT.EQ.0) THEN - B(4)=-B(4) - ELSEIF (IOPT.NE.1) THEN - CALL KTWARN('KTLBST',200,*999) - STOP - ENDIF - DO 30 I=1,4 - DO 20 J=1,4 - C(I,J)=B(I)*B(J)*M - IF (I.EQ.J) C(I,J)=C(I,J)+1 - 20 CONTINUE - 30 CONTINUE - C(4,4)=C(4,4)-2 - CALL KTMMUL(C,R,R) - RETURN - 999 RETURN 1 - END -C----------------------------------------------------------------------- - SUBROUTINE KTRROT(R,A,B,*) - IMPLICIT NONE -C PREMULTIPLY R BY THE 4 BY 4 MATRIX TO -C ROTATE FROM VECTOR A TO VECTOR B BY THE SHORTEST ROUTE -C IF THEY ARE EXACTLY BACK-TO-BACK, THE ROTATION AXIS IS THE VECTOR -C WHICH IS PERPENDICULAR TO THEM AND THE X AXIS, UNLESS THEY ARE -C PERPENDICULAR TO THE Y AXIS, WHEN IT IS THE VECTOR WHICH IS -C PERPENDICULAR TO THEM AND THE Y AXIS. -C NOTE THAT THESE CONDITIONS GUARANTEE THAT IF BOTH ARE PERPENDICULAR -C TO THE Z AXIS, IT WILL BE USED AS THE ROTATION AXIS. -C -C LAST ARGUMENT IS LABEL TO JUMP TO IF EITHER HAS LENGTH ZERO -C - DOUBLE PRECISION R(4,4),M(4,4),A(4),B(4),C(4),D(4),AL,BL,CL,DL,EPS -C---SQRT(2*EPS) IS THE ANGLE IN RADIANS OF THE SMALLEST ALLOWED ROTATION -C NOTE THAT IF YOU CONVERT THIS PROGRAM TO SINGLE PRECISION, YOU WILL -C NEED TO INCREASE EPS TO AROUND 0.5E-4 - PARAMETER (EPS=0.5D-6) - AL=A(1)**2+A(2)**2+A(3)**2 - BL=B(1)**2+B(2)**2+B(3)**2 - IF (AL.LE.0.OR.BL.LE.0) CALL KTWARN('KTRROT',100,*999) - AL=1/SQRT(AL) - BL=1/SQRT(BL) - CL=(A(1)*B(1)+A(2)*B(2)+A(3)*B(3))*AL*BL -C---IF THEY ARE COLLINEAR, DON'T NEED TO DO ANYTHING - IF (CL.GE.1-EPS) THEN - RETURN -C---IF THEY ARE BACK-TO-BACK, USE THE AXIS PERP TO THEM AND X AXIS - ELSEIF (CL.LE.-1+EPS) THEN - IF (ABS(B(2)).GT.EPS) THEN - C(1)= 0 - C(2)=-B(3) - C(3)= B(2) -C---UNLESS THEY ARE PERPENDICULAR TO THE Y AXIS, - ELSE - C(1)= B(3) - C(2)= 0 - C(3)=-B(1) - ENDIF -C---OTHERWISE FIND ROTATION AXIS - ELSE - C(1)=A(2)*B(3)-A(3)*B(2) - C(2)=A(3)*B(1)-A(1)*B(3) - C(3)=A(1)*B(2)-A(2)*B(1) - ENDIF - CL=C(1)**2+C(2)**2+C(3)**2 - IF (CL.LE.0) CALL KTWARN('KTRROT',101,*999) - CL=1/SQRT(CL) -C---FIND ROTATION TO INTERMEDIATE AXES FROM A - D(1)=A(2)*C(3)-A(3)*C(2) - D(2)=A(3)*C(1)-A(1)*C(3) - D(3)=A(1)*C(2)-A(2)*C(1) - DL=AL*CL - M(1,1)=A(1)*AL - M(1,2)=A(2)*AL - M(1,3)=A(3)*AL - M(1,4)=0 - M(2,1)=C(1)*CL - M(2,2)=C(2)*CL - M(2,3)=C(3)*CL - M(2,4)=0 - M(3,1)=D(1)*DL - M(3,2)=D(2)*DL - M(3,3)=D(3)*DL - M(3,4)=0 - M(4,1)=0 - M(4,2)=0 - M(4,3)=0 - M(4,4)=1 - CALL KTMMUL(M,R,R) -C---AND ROTATION FROM INTERMEDIATE AXES TO B - D(1)=B(2)*C(3)-B(3)*C(2) - D(2)=B(3)*C(1)-B(1)*C(3) - D(3)=B(1)*C(2)-B(2)*C(1) - DL=BL*CL - M(1,1)=B(1)*BL - M(2,1)=B(2)*BL - M(3,1)=B(3)*BL - M(1,2)=C(1)*CL - M(2,2)=C(2)*CL - M(3,2)=C(3)*CL - M(1,3)=D(1)*DL - M(2,3)=D(2)*DL - M(3,3)=D(3)*DL - CALL KTMMUL(M,R,R) - RETURN - 999 RETURN 1 - END -C----------------------------------------------------------------------- - SUBROUTINE KTVMUL(M,A,B) - IMPLICIT NONE -C 4 BY 4 MATRIX TIMES 4 VECTOR: B=M*A. -C ALL ARE DOUBLE PRECISION -C IT IS SAFE TO CALL WITH B=A -C FIRST SUBSCRIPT=ROWS, SECOND=COLUMNS - DOUBLE PRECISION M(4,4),A(4),B(4),C(4) - INTEGER I,J - DO 20 I=1,4 - C(I)=0 - DO 10 J=1,4 - C(I)=C(I)+M(I,J)*A(J) - 10 CONTINUE - 20 CONTINUE - DO 30 I=1,4 - B(I)=C(I) - 30 CONTINUE - END -C----------------------------------------------------------------------- - SUBROUTINE KTMMUL(A,B,C) - IMPLICIT NONE -C 4 BY 4 MATRIX MULTIPLICATION: C=A*B. -C ALL ARE DOUBLE PRECISION -C IT IS SAFE TO CALL WITH C=A OR B. -C FIRST SUBSCRIPT=ROWS, SECOND=COLUMNS - DOUBLE PRECISION A(4,4),B(4,4),C(4,4),D(4,4) - INTEGER I,J,K - DO 30 I=1,4 - DO 20 J=1,4 - D(I,J)=0 - DO 10 K=1,4 - D(I,J)=D(I,J)+A(I,K)*B(K,J) - 10 CONTINUE - 20 CONTINUE - 30 CONTINUE - DO 50 I=1,4 - DO 40 J=1,4 - C(I,J)=D(I,J) - 40 CONTINUE - 50 CONTINUE - END -C----------------------------------------------------------------------- - SUBROUTINE KTINVT(A,B) - IMPLICIT NONE -C---INVERT TRANSFORMATION MATRIX A -C -C A = INPUT : 4 BY 4 TRANSFORMATION MATRIX -C B = OUTPUT : INVERTED TRANSFORMATION MATRIX -C -C IF A IS NOT A TRANSFORMATION MATRIX YOU WILL GET STRANGE RESULTS -C -C NOTE THAT IT IS SAFE TO CALL WITH A=B -C - DOUBLE PRECISION A(4,4),B(4,4),C(4,4) - INTEGER I,J -C---TRANSPOSE - DO 20 I=1,4 - DO 10 J=1,4 - C(I,J)=A(J,I) - 10 CONTINUE - 20 CONTINUE -C---NEGATE ENERGY-MOMENTUM MIXING TERMS - DO 30 I=1,3 - C(4,I)=-C(4,I) - C(I,4)=-C(I,4) - 30 CONTINUE -C---OUTPUT - DO 50 I=1,4 - DO 40 J=1,4 - B(I,J)=C(I,J) - 40 CONTINUE - 50 CONTINUE - END -C----------------------------------------------------------------------- - FUNCTION KTMDPI(PHI) - IMPLICIT NONE -C---RETURNS PHI, MOVED ONTO THE RANGE [-PI,PI) - DOUBLE PRECISION KTMDPI,PHI,PI,TWOPI,THRPI,EPS - PARAMETER (PI=3.14159265358979324D0,TWOPI=6.28318530717958648D0, - & THRPI=9.42477796076937972D0) - PARAMETER (EPS=1D-15) - KTMDPI=PHI - IF (KTMDPI.LE.PI) THEN - IF (KTMDPI.GT.-PI) THEN - GOTO 100 - ELSEIF (KTMDPI.GT.-THRPI) THEN - KTMDPI=KTMDPI+TWOPI - ELSE - KTMDPI=-MOD(PI-KTMDPI,TWOPI)+PI - ENDIF - ELSEIF (KTMDPI.LE.THRPI) THEN - KTMDPI=KTMDPI-TWOPI - ELSE - KTMDPI=MOD(PI+KTMDPI,TWOPI)-PI - ENDIF - 100 IF (ABS(KTMDPI).LT.EPS) KTMDPI=0 - END -C----------------------------------------------------------------------- - SUBROUTINE KTWARN(SUBRTN,ICODE,*) -C DEALS WITH ERRORS DURING EXECUTION -C SUBRTN = NAME OF CALLING SUBROUTINE -C ICODE = ERROR CODE: - 99 PRINT WARNING & CONTINUE -C 100-199 PRINT WARNING & JUMP -C 200- PRINT WARNING & STOP DEAD -C----------------------------------------------------------------------- - INTEGER ICODE - CHARACTER*6 SUBRTN - WRITE (6,10) SUBRTN,ICODE - 10 FORMAT(/' KTWARN CALLED FROM SUBPROGRAM ',A6,': CODE =',I4/) - IF (ICODE.LT.100) RETURN - IF (ICODE.LT.200) RETURN 1 - STOP - END -C----------------------------------------------------------------------- -C----------------------------------------------------------------------- -C----------------------------------------------------------------------- diff --git a/GeneratorInterface/LHEInterface/plugins/pathhelper.F b/GeneratorInterface/LHEInterface/plugins/pathhelper.F deleted file mode 100644 index 6a540f1f99de3..0000000000000 --- a/GeneratorInterface/LHEInterface/plugins/pathhelper.F +++ /dev/null @@ -1,31 +0,0 @@ - subroutine mysetPDFpath(pathname) - implicit real*8 (A-H,O-Z) - integer noemax,nopmax,npfmax,nofmax,linemax,nmxset - parameter (noemax=1000,nopmax=40,npfmax=10,nofmax=10,linemax=20) - parameter (nmxset=3) - character*232 pathname - character*232 lhapath - common/LHAPDFC/lhapath - character*20 lhaparm(20) - real*8 lhavalue(20) - common/LHACONTROL/lhaparm,lhavalue - - lhaparm(20) = 'LHAPATH' - do j=1,lnblnk(lhapath) - lhapath(j:j)='' - enddo - lhapath = pathname - - return - end - - subroutine setherwpdf() - CHARACTER*20 AUTPDF - CHARACTER*4 BDECAY - COMMON/HWPRCH/AUTPDF(2),BDECAY - - AUTPDF(1) = 'HWLHAPDF' - AUTPDF(2) = 'HWLHAPDF' - - return - end