Pythia6Interface/plugins/pythia_gluinohadron.f

0001 C*********************************************************************
0002
0003 C...PYRHAD
0004 C...Initialize R-hadron data, like names (for event listing) and
0005 C...charges. So far nothing but the minimum required.
0006 C...Feel free to choose between hadron-like names or
0007 C...flavour-content-based ones.
0008
0009       SUBROUTINE PYGLRHAD
0010
0011 C...Double precision and integer declarations.
0012       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0013 C...Parameter statement to help give large particle numbers
0014 C...(left- and righthanded SUSY, excited fermions).
0015       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KEXCIT=4000000)
0016 C...Commonblocks.
0017       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0018       COMMON/PYDAT4/CHAF(500,2)
0019       CHARACTER CHAF*16
0020       SAVE /PYDAT2/,/PYDAT4/
0021
0022 C...Local R-hadron data arrays, to fill into the global ones.
0023       DIMENSION KFRH(20),KCHGRH(20),KANTRH(20)
0024       CHARACTER*16 CHRHA(20),CHRHB(20)
0025 C...Codes.
0026       DATA KFRH/1000993,1009213,1009313,1009323,1009113,1009223,
0027      &1009333,1091114,1092114,1092214,1092224,1093114,1093214,
0028      &1093224,1093314,1093324,1093334,3*0/
0029 C...Three times charge.
0030       DATA KCHGRH/0,3,0,3,0,0,0,-3,0,3,6,-3,0,3,-3,0,-3,3*0/
0031 C...Existence of a distinct antiparticle.
0032       DATA KANTRH/0,3*1,3*0,10*1,3*0/
0033 C...One possibility: hadron-like names.
0034       DATA CHRHA/'~g glueball','~g rho+','~g K*0','~g K*+',
0035      &'~g rho0','~g omega','~g phi','~g Delta-','~g Delta0',
0036      &'~g Delta+','~g Delta++','~g Sigma*-','~g Sigma*0',
0037      &'~g Sigma*+','~g Xi*-','~g Xi*0 ','~g Omega-',3*' '/
0038       DATA CHRHB/' ','~g rho-','~g K*bar0','~g K*-',3*' ',
0039      &'~g Deltabar+','~g Deltabar0','~g Deltabar-','~g Deltabar--',
0040      &'~g Sigma*bar+','~g Sigma*bar0','~g Sigma*bar-','~g Xi*bar+',
0041      &'~g Xi*bar0','~g Omegabar+',3*' '/
0042 C...Another possibility: flavour-contents-based names.
0043 c      DATA CHRHA/'~g g','~g u dbar','~g d sbar','~g u sbar',
0044 c     &'~g d dbar','~g u ubar','~g s sbar','~g ddd','~g udd',
0045 c     &'~g uud','~g uuu','~g sdd','~g sud','~g suu','~g ssd',
0046 c     &'~g ssu','~g sss',3*' '/
0047 c      DATA CHRHB/' ','~g d ubar','~g s dbar','~g s ubar',3*' ',
0048 c     &'~g ddd bar','~g udd bar','~g uud bar','~g uuu bar',
0049 c     &'~g sdd bar','~g sud bar','~g suu bar','~g ssd bar',
0050 c     &'~g ssu bar','~g sss bar',3*' '/
0051
0052 C...Fill in data.
0053       DO 100 I=1,17
0054         KC=400+I
0055         KCHG(KC,1)=KCHGRH(I)
0056         KCHG(KC,2)=0
0057         KCHG(KC,3)=KANTRH(I)
0058         KCHG(KC,4)=KFRH(I)
0059         CHAF(KC,1)=CHRHA(I)
0060         CHAF(KC,2)=CHRHB(I)
0061   100 CONTINUE
0062
0063       RETURN
0064       END
0065
0066
0067 C*********************************************************************
0068
0069 C...PYGLFR
0070 C...Fragments the string near to a gluino, to form a gluino-hadron,
0071 C...either by producing a new g-g pair or two new q-qbar ones.
0072
0073       SUBROUTINE PYGLFR
0074
0075 C...Double precision and integer declarations.
0076       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0077 C...Parameter statement to help give large particle numbers
0078 C...(left- and righthanded SUSY, excited fermions).
0079       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KEXCIT=4000000)
0080 C...Commonblocks.
0081       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0082       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0083       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0084 C...Note that dimensions below grew from 4000 to 8000 in Pythia 6.2!
0085       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0086       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0087       COMMON/PYINT1/MINT(400),VINT(400)
0088       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0089       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYDAT3/,/PYPARS/,/PYINT1/,
0090      &/PYINT2/
0091 C...Local array.
0092       DIMENSION PSUM(5),KFSAV(2),PMSAV(2),PSAV(2,5)
0093
0094 C...Initialise some variables that may be used before being set
0095       PMNEW=0.0D0
0096       IRECSV=0
0097       IS=0
0098
0099 C...Free parameter: relative probability for gluino-gluon-ball.
0100 C...(But occasional low-mass string will never become it anyway.)
0101       PROBGG=0.1D0
0102
0103 C...Free parameter: gluon constituent mass.
0104       PMGLU=0.7D0
0105
0106 C...Free parameter: max kinetic energy in gluino-hadron.
0107       PMKIN=0.5D0
0108
0109 C...Switch off popcorn baryon production. (Not imperative, but more
0110 C...failed events when popcorn is allowed.)
0111       MSTJ12=MSTJ(12)
0112       MSTJ(12)=1
0113
0114 C...Convenient shorthand.
0115       KFGL=KSUSY1+21
0116
0117 C...Loopback point for serious problems, with new try.
0118       LOOP=0
0119       CALL PYEDIT(21)
0120       CHGSAV=PYP(0,6)
0121    90 LOOP=LOOP+1
0122       IF(LOOP.GT.1) CALL PYEDIT(22)
0123
0124 C...Take copy of string system(s), leaving extra free slot after gluino.
0125 C...(Eventually to be overwritten by one q and one qbar string break.)
0126       NOLD=N
0127       NGLUI=0
0128       DO 120 I=1,NOLD
0129         ICOPY=0
0130         IF(K(I,1).EQ.2) ICOPY=1
0131         IF(K(I,1).EQ.1.AND.I.GE.2) THEN
0132           IF(K(I-1,1).EQ.12) ICOPY=1
0133         ENDIF
0134         IF(ICOPY.EQ.1) THEN
0135           N=N+1
0136           DO 100 J=1,5
0137             K(N,J)=K(I,J)
0138             P(N,J)=P(I,J)
0139             V(N,J)=V(I,J)
0140   100     CONTINUE
0141           K(I,1)=K(I,1)+10
0142           K(I,4)=N
0143           K(I,5)=N
0144           K(N,3)=I
0145           IF(K(I,2).EQ.KFGL) THEN
0146             NGLUI=NGLUI+1
0147             N=N+1
0148             DO 110 J=1,5
0149               K(N,J)=K(N-1,J)
0150               P(N,J)=0D0
0151               V(N,J)=V(I,J)
0152   110       CONTINUE
0153             K(I,5)=N
0154             K(N,2)=21
0155           ENDIF
0156         ENDIF
0157   120 CONTINUE
0158
0159 C...Loop over (up to) two gluinos per event.
0160       DO 300 IGLUI=1,NGLUI
0161
0162 C...Identify position of gluino (randomize order of treatment).
0163         IGL=0
0164         NGL=0
0165         DO 130 I=1,N
0166           IF(K(I,1).EQ.2.AND.K(I,2).EQ.KFGL) THEN
0167             NGL=NGL+1
0168             IF(IGLUI.EQ.NGLUI) THEN
0169               IGL=I
0170             ELSEIF(NGL.EQ.1) THEN
0171               IF(PYR(0).LT.0.5D0) IGL=I
0172             ELSEIF(IGL.EQ.0) THEN
0173               IGL=I
0174             ENDIF
0175           ENDIF
0176   130   CONTINUE
0177
0178 C...Identify range of partons on string the gluino belongs to.
0179         IMIN=IGL
0180   140   IMIN=IMIN-1
0181         IF(K(IMIN-1,1).EQ.2) GOTO 140
0182         IMAX=IGL
0183   150   IMAX=IMAX+1
0184         IF(K(IMAX,1).EQ.2) GOTO 150
0185
0186 C...Find mass of this gluino-string.
0187         DO 170 J=1,5
0188           PSUM(J)=0D0
0189           DO 160 I=IMIN,IMAX
0190             PSUM(J)=PSUM(J)+P(I,J)
0191   160     CONTINUE
0192   170   CONTINUE
0193         PSUM(5)=SQRT(MAX(0D0,PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-
0194      &  PSUM(3)**2))
0195
0196 C...If low-mass, then consider gluino-hadron already formed.
0197         IF(PSUM(5).LE.P(IGL,5)+P(IMIN,5)+P(IMAX,5)+PMKIN) THEN
0198           DO 180 I=IMIN,IMAX
0199             K(I,1)=15+IGLUI
0200   180     CONTINUE
0201           GOTO 300
0202         ENDIF
0203
0204 C...Else break string by production of new gg or two new qqbar pairs.
0205 C...(Also diquarks allowed, but not popcorn, and not two adjacent.)
0206         IF(PYR(0).LT.PROBGG) THEN
0207 C...Let a gluon occupy two slots, to make administration work the same
0208 C...way as for the qqbar case.
0209           KFSAV(1)=21
0210           KFSAV(2)=21
0211           PMSAV(1)=0.5D0*PMGLU
0212           PMSAV(2)=0.5D0*PMGLU
0213         ELSE
0214   185     CALL PYDCYK(K(IMIN,2),0,KFSAV(1),KFDUM)
0215           CALL PYDCYK(K(IMAX,2),0,KFSAV(2),KFDUM)
0216           IF(IABS(KFSAV(1)).GT.10.AND.IABS(KFSAV(2)).GT.10) GOTO 185
0217           IF(IABS(KFSAV(1)).GT.10.AND.IABS(K(IGL-1,2)).GT.10) GOTO 185
0218           IF(IABS(KFSAV(2)).GT.10.AND.IABS(K(IGL+2,2)).GT.10) GOTO 185
0219           KFSAV(1)=ISIGN(MOD(IABS(KFSAV(1)),10000),KFSAV(1))
0220           KFSAV(2)=ISIGN(MOD(IABS(KFSAV(2)),10000),KFSAV(2))
0221           MSTJ(93)=1
0222           PMSAV(1)=PYMASS(KFSAV(1))
0223           MSTJ(93)=1
0224           PMSAV(2)=PYMASS(KFSAV(2))
0225         ENDIF
0226
0227 C...Mass of gluino-hadron.
0228         PMGSAV=P(IGL,5)
0229         PMGB=P(IGL,5)+PMSAV(1)+PMSAV(2)
0230
0231 C...Pick at random order in which both sides of gluino string break.
0232         ISIDE=INT(1.5D0+PYR(0))
0233         DO 220 ISDE=1,2
0234           IF(ISDE.EQ.1) IS=ISIDE
0235           IF(ISDE.EQ.2) IS=3-ISIDE
0236
0237 C...Pick momentum sharing according to fragmentation function as if bottom.
0238           PMBSAV=PARF(105)
0239           PARF(105)=PMGSAV
0240           CALL PYZDIS(5,0,PMGB**2,ZGL)
0241           PARF(105)=PMBSAV
0242           ZGL=MAX(0.9D0,MIN(0.9999D0,ZGL))
0243           DO 190 J=1,5
0244             PSAV(IS,J)=(1D0-ZGL)*P(IGL,J)
0245             P(IGL,J)=ZGL*P(IGL,J)
0246   190    CONTINUE
0247
0248 C...Target gluino-hadron mass for this stage of momentum reshuffling.
0249           PMOLD=P(IGL,5)
0250           IF(ISDE.EQ.1) PMNEW=PMGSAV+PMSAV(IS)
0251           IF(ISDE.EQ.2) PMNEW=PMGB
0252
0253 C...Recoiling parton from which to shuffle momentum. System momentum.
0254           IF(IS.EQ.1) IREC=IGL-1
0255           IF(IS.EQ.2) IREC=IGL+2
0256   200     DO 210 J=1,4
0257             PSUM(J)=P(IGL,J)+P(IREC,J)
0258   210     CONTINUE
0259
0260 C...Boost to rest frame of system, and align gluino along +z axis.
0261           CALL PYROBO(IGL,IGL,0D0,0D0,-PSUM(1)/PSUM(4),
0262      &    -PSUM(2)/PSUM(4),-PSUM(3)/PSUM(4))
0263           CALL PYROBO(IREC,IREC,0D0,0D0,-PSUM(1)/PSUM(4),
0264      &    -PSUM(2)/PSUM(4),-PSUM(3)/PSUM(4))
0265           PHI=PYANGL(P(IGL,1),P(IGL,2))
0266           CALL PYROBO(IGL,IGL,0D0,-PHI,0D0,0D0,0D0)
0267           CALL PYROBO(IREC,IREC,0D0,-PHI,0D0,0D0,0D0)
0268           THETA=PYANGL(P(IGL,3),P(IGL,1))
0269           CALL PYROBO(IGL,IGL,-THETA,0D0,0D0,0D0,0D0)
0270           CALL PYROBO(IREC,IREC,-THETA,0D0,0D0,0D0,0D0)
0271
0272 C...Calculate new kinematics in this frame, for desired gluino mass.
0273           ETOT=P(IGL,4)+P(IREC,4)
0274           IF(ETOT.GT.PMNEW+P(IREC,5)) THEN
0275             IFAIL=0
0276             PZNEW=0.5D0*SQRT(MAX(0D0,(ETOT**2-PMNEW**2-P(IREC,5)**2)**2-
0277      &      4D0*PMNEW**2*P(IREC,5)**2))/ETOT
0278             P(IGL,3)=PZNEW
0279             P(IGL,4)=SQRT(PZNEW**2+PMNEW**2)
0280             P(IGL,5)=PMNEW
0281             P(IREC,3)=-PZNEW
0282             P(IREC,4)=SQRT(PZNEW**2+P(IREC,5)**2)
0283
0284 C...If not enough momentum, take what can be taken.
0285           ELSE
0286             IFAIL=1
0287             PMOLD=ETOT-P(IREC,5)
0288             P(IGL,3)=0D0
0289             P(IGL,4)=PMOLD
0290             P(IGL,5)=PMOLD
0291             P(IREC,3)=0D0
0292             P(IREC,4)=P(IREC,5)
0293           ENDIF
0294
0295 C...Bost back to lab frame.
0296           CALL PYROBO(IGL,IGL,THETA,PHI,PSUM(1)/PSUM(4),
0297      &    PSUM(2)/PSUM(4),PSUM(3)/PSUM(4))
0298           CALL PYROBO(IREC,IREC,THETA,PHI,PSUM(1)/PSUM(4),
0299      &    PSUM(2)/PSUM(4),PSUM(3)/PSUM(4))
0300
0301 C...Loop back when not enough momentum could be shuffled.
0302 C...(As long as there is something left on either side.)
0303           IF(IFAIL.EQ.1) THEN
0304   215       IF(IS.EQ.1.AND.IREC.GT.IMIN) THEN
0305               IREC=IREC-1
0306               GOTO 200
0307             ELSEIF(IS.EQ.2.AND.IREC.LT.IMAX) THEN
0308               IREC=IREC+1
0309               GOTO 200
0310             ELSEIF(ISDE.EQ.2.AND.IS.EQ.3-ISIDE) THEN
0311               IS=ISIDE
0312               IREC=IRECSV
0313               GOTO 215
0314             ENDIF
0315           ENDIF
0316
0317 C...End loop over fragmentation of two sides around gluino.
0318          IRECSV=IREC
0319   220   CONTINUE
0320
0321 C...New slot at end of record for gluino R-hadron.
0322         DO 230 J=1,5
0323           K(N+1,J)=0
0324           P(N+1,J)=P(IGL,J)
0325           V(N+1,J)=V(IGL,J)
0326   230   CONTINUE
0327
0328 C...Status and code of this slot.
0329         K(N+1,1)=15+IGLUI
0330         KFSVMX=MAX(IABS(KFSAV(1)),IABS(KFSAV(2)))
0331         KFSVMN=MIN(IABS(KFSAV(1)),IABS(KFSAV(2)))
0332 C...Gluino-ball.
0333         IF(KFSVMX.EQ.21) THEN
0334           K(N+1,2)=KSUSY1+993
0335 C...Gluino-meson.
0336         ELSEIF(KFSVMX.LT.10) THEN
0337           K(N+1,2)=KSUSY1+9000+100*KFSVMX+10*KFSVMN+3
0338           IF(KFSVMX.EQ.KFSVMN) THEN
0339           ELSEIF(MOD(KFSVMX,2).EQ.0) THEN
0340             IF(KFSVMX.EQ.KFSAV(1).OR.KFSVMX.EQ.KFSAV(2))
0341      &      K(N+1,2)=-K(N+1,2)
0342           ELSE
0343             IF(KFSVMX.EQ.-KFSAV(1).OR.KFSVMX.EQ.-KFSAV(2))
0344      &      K(N+1,2)=-K(N+1,2)
0345           ENDIF
0346 C...Gluino-baryon.
0347         ELSE
0348           KFSVX1=KFSVMX/1000
0349           KFSVX2=MOD(KFSVMX/100,10)
0350           KFA=MAX(KFSVX1,KFSVX2,KFSVMN)
0351           KFC=MIN(KFSVX1,KFSVX2,KFSVMN)
0352           KFB=KFSVX1+KFSVX2+KFSVMN-KFA-KFC
0353           K(N+1,2)=SIGN(KSUSY1+90000+1000*KFA+100*KFB+10*KFC+4,
0354      &    -KFSAV(1))
0355         ENDIF
0356         K(N+1,3)=K(IGL,3)
0357         N=N+1
0358
0359 C...Code and momentum of two new string endpoints.
0360         K(IGL,2)=KFSAV(1)
0361         K(IGL+1,2)=KFSAV(2)
0362         IF(KFSAV(1).NE.21) K(IGL,1)=1
0363         DO 240 J=1,5
0364           P(IGL,J)=PSAV(1,J)
0365           P(IGL+1,J)=PSAV(2,J)
0366   240   CONTINUE
0367
0368 C...End of loop over two gluinos.
0369   300 CONTINUE
0370
0371 C...Cleanup: remove zero-energy gluons.
0372       NNOW=N
0373       N=NOLD
0374       DO 330 I=NOLD+1,NNOW
0375         IF(K(I,2).EQ.21.AND.P(I,4).LT.1D-10) THEN
0376         ELSEIF(I.EQ.N+1) THEN
0377           N=N+1
0378         ELSE
0379           N=N+1
0380           DO 320 J=1,5
0381             K(N,J)=K(I,J)
0382             P(N,J)=P(I,J)
0383             V(N,J)=V(I,J)
0384   320     CONTINUE
0385         ENDIF
0386   330 CONTINUE
0387
0388 C...Finish off with standard hadronization.
0389 C...Note that the R-hadrons are not affected here, since they
0390 C...(deliberately erroneously) have been bookkept as decayed.
0391       CALL PYEXEC
0392
0393 C...Restore code of gluino-hadrons to non-fragmented numbers.
0394       N6=0
0395       N7=0
0396       DO 340 I=1,N
0397         IF(K(I,1).EQ.16.OR.K(I,1).EQ.17) K(I,1)=K(I,1)-10
0398         IF(K(I,1).EQ.6) N6=N6+1
0399         IF(K(I,1).EQ.7) N7=N7+1
0400   340 CONTINUE
0401       IF(N6.GT.1.OR.N7.GT.1) MSTU(24)=1
0402
0403 C...Extracheck charge.
0404       CHGNEW=PYP(0,6)
0405       IF(ABS(CHGNEW-CHGSAV).GT.0.1D0) MSTU(24)=1
0406
0407 C...In case of trouble, make up to five attempts.
0408       IF(MSTU(24).NE.0.AND.LOOP.LT.5) THEN
0409         WRITE(*,*) '     ...give it new try...'
0410         MSTU(23)=MSTU(23)-1
0411         GOTO 90
0412       ELSEIF(MSTU(24).NE.0) THEN
0413         WRITE(*,*) '     ...but still fail after repeated attempts!'
0414       ELSEIF(MSTU(24).EQ.0.AND.LOOP.GT.1) THEN
0415         WRITE(*,*) '     ...and now it worked!'
0416       ENDIF
0417
0418 C...Finished. Restore baryon production model.
0419       MSTJ(12)=MSTJ12
0420
0421       RETURN
0422       END
0423