Pythia6Interface/plugins/pythia_stophadron.f

0001 C*********************************************************************
0002
0003 C...PYRHAD
0004 C...Initialize stop_1 R-hadron data, like names (for event listing)
0005 C...and 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 PYSTRHAD
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/1000612,1000622,1000632,1000642,1000652,1006113,
0027      &1006211,1006213,1006223,1006311,1006313,1006321,1006323,
0028      &1006333,6*0/
0029 C...Three times charge.
0030       DATA KCHGRH/3,0,3,0,3,0,3,3,6,0,0,3,3,0,6*0/
0031 C...Existence of a distinct antiparticle.
0032       DATA KANTRH/14*1,6*0/
0033 C...One possibility: hadron-like names.
0034       DATA CHRHA/'~T+','~T0','~T_s+','~T_c0','~T_b+','~T_dd10',
0035      &'~T_ud0+','~T_ud1+','~T_uu1++','~T_sd00','~T_sd10',
0036      &'~T_su0+','~T_su1+','~T_ss10',6*' '/
0037       DATA CHRHB/'~Tbar-','~Tbar0','~Tbar_s-','~Tbar_c0','~Tbar_b-',
0038      &'~Tbar_dd10','~Tbar_ud0-','~Tbar_ud1-','~Tbar_uu1--',
0039      &'~Tbar_sd00','~Tbar_sd10','~Tbar_su0-','~Tbar_su1-',
0040      &'~Tbar_ss10',6*' '/
0041 C...Another possibility: flavour-contents-based names.
0042 C      DATA CHRHA/'~t dbar','~t ubar','~t sbar','~t cbar','~t bbar',
0043 C     &'~t dd1','~t ud0','~t ud1','~t uu1','~t sd0','~t sd1','~t su0',
0044 C     &'~t su1','~t ss1',6*' '/
0045 C      DATA CHRHB/'~tbar d','~tbar u','~tbar s','~tbar c','~tbar b',
0046 C     &'~tbar dd1bar','~tbar ud0bar','~tbar ud1bar','~tbar uu1bar',
0047 C     &'~tbar sd0bar','~tbar sd1bar','~tbar su0bar','~tbar su1bar',
0048 C     &'~tbar ss1bar',6*' '/
0049
0050 C...Fill in data.
0051       DO 100 I=1,14
0052         KC=400+I
0053         KCHG(KC,1)=KCHGRH(I)
0054         KCHG(KC,2)=0
0055         KCHG(KC,3)=KANTRH(I)
0056         KCHG(KC,4)=KFRH(I)
0057         CHAF(KC,1)=CHRHA(I)
0058         CHAF(KC,2)=CHRHB(I)
0059   100 CONTINUE
0060
0061       RETURN
0062       END
0063
0064 C*********************************************************************
0065
0066 C...PYSTFR
0067 C...Fragments the string near to a stop, to form a stop-hadron,
0068 C...by producing a new q-qbar pair.
0069
0070       SUBROUTINE PYSTFR(IERR)
0071
0072 C...Double precision and integer declarations.
0073       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0074       INTEGER PYCOMP
0075 C...Parameter statement to help give large particle numbers
0076 C...(left- and righthanded SUSY, excited fermions).
0077       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KEXCIT=4000000)
0078 C...Commonblocks.
0079       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0080       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0081       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0082 C...Note that dimensions below grew from 4000 to 8000 in Pythia 6.2!
0083       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0084       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0085       COMMON/PYINT1/MINT(400),VINT(400)
0086       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0087       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYDAT3/,/PYPARS/,/PYINT1/,
0088      &/PYINT2/
0089 C...Local array.
0090       DIMENSION PSUM(5),PSAV(5),IPOSST(10)
0091
0092 C...Default is no error.
0093       IERR=0
0094
0095 C...Free parameter: max kinetic energy in gluino-hadron.
0096       PMKIN=0.5D0
0097
0098 C...Free parameter: part of stop mass that does not participate
0099 C...in weak decay.
0100       PMINAC=0.5D0
0101
0102 C...Switch off popcorn baryon production. (Not imperative, but more
0103 C...failed events when popcorn is allowed.)
0104       MSTJ12=MSTJ(12)
0105       MSTJ(12)=1
0106
0107 C...Convenient shorthand.
0108       KFST=KSUSY1+6
0109       KCST=PYCOMP(KFST)
0110       KFGL=KSUSY1+21
0111
0112 C...Loopback point for serious problems, with new try.
0113       LOOP=0
0114       CALL PYEDIT(21)
0115       CHGSAV=PYP(0,6)
0116    90 LOOP=LOOP+1
0117       IF(LOOP.GT.1) CALL PYEDIT(22)
0118
0119 C...Give up when too much problems.
0120       IF(LOOP.GT.5) THEN
0121         WRITE(*,*) ' Problematical event skipped'
0122         IERR=1
0123         RETURN
0124       ENDIF
0125
0126 C...Take copy of string system(s).
0127       NOLD=N
0128       NSTOP=0
0129       DO 120 I=1,NOLD
0130         ICOPY=0
0131         IF(K(I,1).EQ.2) ICOPY=1
0132         IF(K(I,1).EQ.1.AND.I.GE.2) THEN
0133           IF(K(I-1,1).EQ.12) ICOPY=1
0134         ENDIF
0135         IF(ICOPY.EQ.1) THEN
0136           N=N+1
0137           DO 100 J=1,5
0138             K(N,J)=K(I,J)
0139             P(N,J)=P(I,J)
0140             V(N,J)=V(I,J)
0141   100     CONTINUE
0142           K(I,1)=K(I,1)+10
0143           K(I,4)=N
0144           K(I,5)=N
0145           K(N,3)=I
0146           IF(IABS(K(I,2)).EQ.KFST) THEN
0147             NSTOP=NSTOP+1
0148             IPOSST(NSTOP)=N
0149           ENDIF
0150         ENDIF
0151   120 CONTINUE
0152       NTMP=N
0153
0154 C...Loop over (up to) two stops per event.
0155 C...Identify position of stop (randomize order of treatment).
0156       IRNST=INT(1.5D0+PYR(0))
0157       DO 300 ISTOP=1,NSTOP
0158         IST=IPOSST(1)
0159         IF(NSTOP.EQ.2.AND.ISTOP.NE.IRNST) IST=IPOSST(2)
0160
0161 C...Identify range of partons on string the stop belongs to.
0162         IMIN=IST+1
0163   140   IMIN=IMIN-1
0164         IF(K(IMIN-1,1).EQ.2) GOTO 140
0165         IMAX=IST-1
0166   150   IMAX=IMAX+1
0167         IF(K(IMAX,1).EQ.2) GOTO 150
0168         IOTHER=IMAX
0169         IF(IST.EQ.IMAX) IOTHER=IMIN
0170
0171 C...Find mass of this stop-string.
0172         DO 170 J=1,5
0173           PSUM(J)=0D0
0174           DO 160 I=IMIN,IMAX
0175             PSUM(J)=PSUM(J)+P(I,J)
0176   160     CONTINUE
0177   170   CONTINUE
0178         PSUM(5)=SQRT(MAX(0D0,PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-
0179      &  PSUM(3)**2))
0180
0181 C...If low-mass, then consider stop-hadron already formed.
0182         IF(PSUM(5).LE.P(IST,5)+P(IOTHER,5)+PMKIN) THEN
0183           DO 180 I=IMIN,IMAX
0184             K(I,1)=2
0185             IF(I.EQ.IMAX) K(I,1)=1
0186             IF(I.NE.IST) THEN
0187               DO 175 J=1,5
0188                 P(IST,J)=P(IST,J)+P(I,J)
0189                 P(I,J)=0D0
0190   175         CONTINUE
0191             ENDIF
0192   180     CONTINUE
0193           P(IST,5)=SQRT(MAX(0D0,P(IST,4)**2-P(IST,1)**2-P(IST,2)**2-
0194      &    P(IST,3)**2))
0195           GOTO 300
0196         ENDIF
0197
0198 C...Else break string by production of new qqbar pair.
0199 C...(Also diquarks allowed, but not popcorn.)
0200         INFLAV=ISIGN(4,K(IST,2))
0201         CALL PYDCYK(INFLAV,0,KFSAV,KFDUM)
0202         KFSAV=ISIGN(MOD(IABS(KFSAV),10000),KFSAV)
0203         MSTJ(93)=1
0204         PMSAV=PYMASS(KFSAV)
0205
0206 C...Mass of stop-hadron.
0207         PMSSAV=P(IST,5)
0208         PMSHAD=P(IST,5)+PMSAV
0209
0210 C...Pick momentum sharing according to fragmentation function as if bottom.
0211         PMBSAV=PARF(105)
0212         PARF(105)=PMSSAV
0213         CALL PYZDIS(5,0,PMSHAD**2,ZST)
0214         PARF(105)=PMBSAV
0215         ZST=MAX(0.9D0,MIN(0.9999D0,ZST))
0216         DO 190 J=1,5
0217           PSAV(J)=(1D0-ZST)*P(IST,J)
0218           P(IST,J)=ZST*P(IST,J)
0219   190  CONTINUE
0220
0221 C...Recoiling parton from which to shuffle momentum. System momentum.
0222         IF(IST.EQ.IMIN) IREC=IST+1
0223         IF(IST.EQ.IMAX) IREC=IST-1
0224   200   DO 210 J=1,4
0225           PSUM(J)=P(IST,J)+P(IREC,J)
0226   210   CONTINUE
0227
0228 C...Boost to rest frame of system, and align stop along +z axis.
0229         CALL PYROBO(IST,IST,0D0,0D0,-PSUM(1)/PSUM(4),
0230      &  -PSUM(2)/PSUM(4),-PSUM(3)/PSUM(4))
0231         CALL PYROBO(IREC,IREC,0D0,0D0,-PSUM(1)/PSUM(4),
0232      &  -PSUM(2)/PSUM(4),-PSUM(3)/PSUM(4))
0233         PHI=PYANGL(P(IST,1),P(IST,2))
0234         CALL PYROBO(IST,IST,0D0,-PHI,0D0,0D0,0D0)
0235         CALL PYROBO(IREC,IREC,0D0,-PHI,0D0,0D0,0D0)
0236         THETA=PYANGL(P(IST,3),P(IST,1))
0237         CALL PYROBO(IST,IST,-THETA,0D0,0D0,0D0,0D0)
0238         CALL PYROBO(IREC,IREC,-THETA,0D0,0D0,0D0,0D0)
0239
0240 C...Calculate new kinematics in this frame, for desired stop hadron mass.
0241         ETOT=P(IST,4)+P(IREC,4)
0242         PMREC=P(IREC,5)
0243         IF(K(IREC,2).NE.21.AND.IABS(K(IREC,2)).NE.KFST) THEN
0244           MSTJ(93)=1
0245           PMREC=PYMASS(K(IREC,2))
0246         ENDIF
0247         IF(ETOT.GT.PMSHAD+PMREC) THEN
0248           IFAIL=0
0249           PZNEW=0.5D0*SQRT(MAX(0D0,(ETOT**2-PMSHAD**2-PMREC**2)**2-
0250      &    4D0*PMSHAD**2*PMREC**2))/ETOT
0251           P(IST,3)=PZNEW
0252           P(IST,4)=SQRT(PZNEW**2+PMSHAD**2)
0253           P(IST,5)=PMSHAD
0254           P(IREC,3)=-PZNEW
0255           P(IREC,4)=SQRT(PZNEW**2+PMREC**2)
0256           P(IREC,5)=PMREC
0257
0258 C...If not enough momentum, take what can be taken.
0259         ELSE
0260           IFAIL=1
0261           P(IST,3)=0D0
0262           P(IST,4)=ETOT-PMREC
0263           P(IST,5)=P(IST,4)
0264           P(IREC,3)=0D0
0265           P(IREC,4)=PMREC
0266           P(IREC,5)=PMREC
0267         ENDIF
0268
0269 C...Bost back to lab frame.
0270         CALL PYROBO(IST,IST,THETA,PHI,PSUM(1)/PSUM(4),
0271      &  PSUM(2)/PSUM(4),PSUM(3)/PSUM(4))
0272         CALL PYROBO(IREC,IREC,THETA,PHI,PSUM(1)/PSUM(4),
0273      &  PSUM(2)/PSUM(4),PSUM(3)/PSUM(4))
0274
0275 C...Loop back when not enough momentum could be shuffled.
0276 C...(As long as there is something left.)
0277         IF(IFAIL.EQ.1) THEN
0278           IF(IST.EQ.IMIN.AND.IREC.LT.IMAX) THEN
0279             IREC=IREC+1
0280             GOTO 200
0281           ELSEIF(IST.EQ.IMAX.AND.IREC.GT.IMIN) THEN
0282             IREC=IREC-1
0283             GOTO 200
0284           ENDIF
0285         ENDIF
0286
0287 C...Particle code for stop-hadron.
0288         KFSTHD=0
0289         IF(K(IST,2).GT.0) THEN
0290           IF(KFSAV.LE.-1.AND.KFSAV.GE.-5) KFSTHD=KSUSY1+600-10*KFSAV+2
0291           IF(KFSAV.GE.1103.AND.KFSAV.LE.3303) KFSTHD=KSUSY1+6000+
0292      &    (KFSAV/10)+MOD(KFSAV,10)
0293         ELSE
0294           IF(KFSAV.GE.1.AND.KFSAV.LE.5) KFSTHD=KSUSY1+600+10*KFSAV+2
0295           IF(KFSAV.LE.-1103.AND.KFSAV.GE.-3303) KFSTHD=KSUSY1+6000+
0296      &    (IABS(KFSAV)/10)+MOD(IABS(KFSAV),10)
0297           KFSTHD=-KFSTHD
0298         ENDIF
0299         IF(KFSTHD.EQ.0) THEN
0300           WRITE(*,*) ' Failed to find R-hadron code from ',
0301      &    K(IST,2),KFSAV
0302           IERR=1
0303           RETURN
0304         ENDIF
0305
0306 C...New slot at end of record for stop-hadron
0307         DO 230 J=1,5
0308           K(N+1,J)=0
0309           P(N+1,J)=P(IST,J)
0310           V(N+1,J)=V(IST,J)
0311   230   CONTINUE
0312         K(N+1,1)=5+ISTOP
0313         K(N+1,2)=KFSTHD
0314         K(N+1,3)=K(IST,3)
0315         N=N+1
0316
0317 C...Code and momentum of new string endpoint.
0318         K(IST,2)=-KFSAV
0319         DO 240 J=1,5
0320           P(IST,J)=PSAV(J)
0321   240   CONTINUE
0322
0323 C...End of loop over two stops.
0324   300 CONTINUE
0325
0326 C...Cleanup: remove zero-energy gluons.
0327       NNOW=N
0328       N=NOLD
0329       DO 330 I=NOLD+1,NNOW
0330         IF(K(I,2).EQ.21.AND.P(I,4).LT.1D-10) THEN
0331         ELSEIF(I.EQ.N+1) THEN
0332           N=N+1
0333         ELSE
0334           N=N+1
0335           DO 320 J=1,5
0336             K(N,J)=K(I,J)
0337             P(N,J)=P(I,J)
0338             V(N,J)=V(I,J)
0339   320     CONTINUE
0340         ENDIF
0341   330 CONTINUE
0342       NNOW=N
0343
0344 C...Check that no low-mass system of diquark-antidiquark kind,
0345 C...or very low-mass of any kind.
0346       KFBEG=0
0347       DO 332 J=1,5
0348         PSUM(J)=0D0
0349   332 CONTINUE
0350       DO 338 I=NOLD+1,NNOW
0351         DO 334 J=1,4
0352           PSUM(J)=PSUM(J)+P(I,J)
0353   334   CONTINUE
0354         IF(KFBEG.EQ.0) THEN
0355           KFBEG=IABS(K(I,2))
0356           MSTJ(93)=1
0357           PSUM(5)=PSUM(5)+PYMASS(K(I,2))
0358         ELSEIF(K(I,1).EQ.1) THEN
0359           KFEND=IABS(K(I,2))
0360           MSTJ(93)=1
0361           PSUM(5)=PSUM(5)+PYMASS(K(I,2))
0362           DELTA=SQRT(MAX(0D0,PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-
0363      &    PSUM(3)**2))-PSUM(5)
0364           IF(KFBEG.GT.10.AND.KFBEG.LT.10000.AND.KFEND.GT.10.AND.
0365      &    KFEND.LT.10000.AND.DELTA.LT.PARJ(32).AND.(KFBEG.NE.21
0366      &    .AND.KFEND.NE.21)) GOTO 90
0367           IF(DELTA.LT.0D0) GOTO 90
0368           KFBEG=0
0369           DO 336 J=1,5
0370             PSUM(J)=0D0
0371   336     CONTINUE
0372         ENDIF
0373   338 CONTINUE
0374
0375 C...Finished with stop hadronization. Restore baryon production model.
0376       MSTJ(12)=MSTJ12
0377
0378 C...Now hadronize everything else. Some cheating to allow sensible
0379 C...momentum shuffling.
0380       MSTJ16=MSTJ(16)
0381       MSTJ(16)=0
0382       CALL PYEXEC
0383       MSTJ(16)=MSTJ16
0384       IF(MSTU(24).NE.0) THEN
0385         WRITE(*,*) ' Event to be skipped'
0386         IERR=1
0387       ENDIF
0388
0389       RETURN
0390       END