* * $Id: decay.F,v 1.1.1.1 1996/01/11 14:14:34 mclareni Exp $ * * $Log: decay.F,v $ * Revision 1.1.1.1 1996/01/11 14:14:34 mclareni * Cojets * * #include "cojets/pilot.h" SUBROUTINE DECAY(IPD,I) C *********************** C-- DECAYS PARTICLES GENERATED BY JETGEN C-- REVISED: 88/04/27 #if defined(CERNLIB_SINGLE) IMPLICIT REAL (A-H,O-Z) #endif #if defined(CERNLIB_DOUBLE) IMPLICIT DOUBLE PRECISION (A-H,O-Z) #endif #include "cojets/data1.inc" #include "cojets/data2.inc" #include "cojets/decpar.inc" #include "cojets/edpar.inc" #include "cojets/event.inc" #include "cojets/forcsl.inc" #include "cojets/itapes.inc" #include "cojets/jet.inc" #include "cojets/maxn.inc" #include "cojets/par.inc" #include "cojets/stable.inc" #include "cojets/zwpar.inc" DIMENSION PM(10,5), RND(10), U(3), BE(3) INTEGER IFLQ(4),IFLR(4) INTEGER IQV(3) PAWT(A,B,C)=SQRT((A**2-(B+C)**2)*(A**2-(B-C)**2))/(2.*A) FOUR(I,J)=P(I,4)*P(J,4)-P(I,1)*P(J,1)-P(I,2)*P(J,2)-P(I,3)*P(J,3) C IF(I+5.GT.MAXJTP) GO TO 500 C 1 DECAY CHANNEL CHOICE, GIVES DECAY PRODUCTS TBR=CJRN(0) KIPD=K(IPD,2) IF(IFORSL.EQ.1.AND.RDECAY(KIPD).LT.0.) THEN C-- FORCE SL IF(CJRN(0.)*BRXFSL.GT.-MOD(RDECAY(KIPD),1.)) THEN C-- REJECT EVENT (CORRECTION WITH RESPECT TO GLOBAL MAX SL BR) IFORSL=-1 RETURN ENDIF IFR=1 IDC=-INT(RDECAY(KIPD))-1 ELSE IFR=0 IDC=IDB(K(IPD,2))-1 ENDIF 100 IDC=IDC+1 IF(IFR.EQ.1) THEN C-- FORCE SL IF(TBR.GT.CBRF(IDC)) GO TO 100 IDC=IKDP(IDC) ELSE IF(TBR.GT.CBR(IDC)) GOTO 100 ENDIF NP=0 IQ=0 PSUMP=0. I1LW=I+1 I1UP=I+5 DO 110 I1=I1LW,I1UP K(I1,2)=KDP(IDC,I1-I) IF(K(I1,2).GT.0.AND.K(I1,2).LT.1000) THEN NP=NP+1 K(I1,1)=-IPD P(I1,5)=PMAS(K(I1,2)) PSUMP=PSUMP+P(I1,5) ELSEIF(K(I1,2).NE.0) THEN C-- Q-QBAR OR DIQUARK-QUARK DECAYS -- STANDARD INTERNAL CODE IQ=IQ+1 IFLQ(IQ)=K(I1,2)/1000 ENDIF 110 CONTINUE NQ=IQ/2 C 2 DECAY MULTIPLICITY ND=NP+NQ PSUM=PSUMP IF(NQ.EQ.0) GOTO 170 119 IF(NP+NQ.EQ.3.AND.K(I+1,2).GE.7.AND.K(I+1,2).LE.18) GO TO 121 CND=CND1*LOG(P(IPD,5)/CND2) 120 GAUSS=SQRT(-2.*CND*LOG(CJRN(1)))*SIN(6.2832*CJRN(2)) ND=0.5+0.5*(NP+NQ)+CND+GAUSS IF(I+ND.GT.MAXJTP) GO TO 500 IF(ND.LT.NP+NQ.OR.ND.LT.2.OR.ND.GT.10) GOTO 120 121 CONTINUE IF(I+ND.GT.MAXJTP) GO TO 500 C-- GENERATION AND PAIRING OFF OF Q-QB AND Q-DIQUARK PAIRS C-- GENERATION AND PAIRING OFF OF Q-QB AND Q-DIQUARK PAIRS DO 130 J=1,4 130 IFLR(J)=IFLQ(J) IF(ND.EQ.NP+NQ) GOTO 150 I1LW=I+NP+1 I1UP=I+ND-NQ DO 140 I1=I1LW,I1UP IQ=1+INT(2*NQ*CJRN(0)) IFLN=SIGN(1+INT(CJRN(0)/PUD),-IFLR(IQ)) IF(ABS(IFLR(IQ)).GT.10) IFLN=-IFLN K(I1,1)=KRCOMB(IFLR(IQ),IFLN) 140 IFLR(IQ)=-IFLN 150 IQ=2+2*INT(NQ*CJRN(0)) K(I+ND-NQ+1,1)=KRCOMB(IFLR(1),IFLR(IQ)) IF(NQ.EQ.2) K(I+ND,1)=KRCOMB(IFLR(3),IFLR(6-IQ)) C-- PARTICLES ARE FORMED, IF TOTAL MASS TOO LARGE ITERATE PSUM=PSUMP I1LW=I+NP+1 I1UP=I+ND DO 160 I1=I1LW,I1UP IF(ABS(K(I1,1)).LT.100) THEN ISPIN=PS1(MOD(ABS(IDENTF(20+K(I1,1))),100)/10)+CJRN(0) K(I1,2)=20+40*ISPIN+K(I1,1) IF(K(I1,1).GE.31.AND.K(I1,1).LE.33) THEN TMIX=CJRN(0) KM=K(I1,1)-30+3*ISPIN K(I1,2)=51+40*ISPIN+INT(TMIX+CMIX(KM,1))+INT(TMIX+CMIX(KM,2)) ENDIF ELSE KA=ABS(K(I1,1)) IQV(1)=KA/100 IQV(2)=MOD(KA/10,10) IQV(3)=MOD(KA,10) CALL BARYO(IQV,ISBARY) K(I1,2)=INTID(SIGN(ISBARY,K(I1,1))) ENDIF K(I1,1)=-IPD P(I1,5)=PMAS(K(I1,2)) 160 PSUM=PSUM+P(I1,5) IF(PSUM+.005.GT.P(IPD,5)) GOTO 119 C 5 PRODUCT CHARGE CONJUGATION FOR DECAYING ANTIPARTICLE C-- BZD,BZS MIXINGS 170 CONTINUE IDEN=IDENTF(K(IPD,2)) IF(ABS(IDEN).EQ.260.AND.CJRN(0.).LT.BZDMIX) IDEN=-IDEN IF(ABS(IDEN).EQ.360.AND.CJRN(0.).LT.BZSMIX) IDEN=-IDEN IF(IDEN.LT.0) THEN I1LW=I+1 I1UP=I+ND DO 180 I1=I1LW,I1UP IDINT=K(I1,2) IDENA=ABS(IDENTF(IDINT)) IF(IDENA.LE.100) THEN IF(LCHARG(IDINT).EQ.0) GO TO 180 ELSEIF(IDENA.LT.1000) THEN IQ1=IDENA/100 IQ2=MOD(IDENA/10,10) IF(IQ1.EQ.IQ2) GO TO 180 ENDIF K(I1,2)=INTID(-IDENTF(IDINT)) 180 CONTINUE ENDIF C 6 ONE-PARTICLE DECAY, START ND-PARTICLE DECAY IF(ND.EQ.1) THEN DO 190 J=1,4 190 P(I+1,J)=P(IPD,J) GOTO 330 ENDIF DO 200 J=1,5 200 PM(1,J)=P(IPD,J) PM(ND,5)=P(I+ND,5) IF(ND.EQ.2) GOTO 260 C-- KROLL-WADA MASS DISTRIBUTION FOR PI0, ETA DALITZ DECAYS IF((K(IPD,2).EQ.51.OR.K(IPD,2).EQ.52).AND. 1(K(I+2,2).EQ.7.OR.K(I+2,2).EQ.8)) THEN E2M=2.*PMAS(7) 201 EEM=E2M*(P(IPD,5)/E2M)**CJRN(0) REEM=(E2M/EEM)**2 WTEE=(1.-(EEM/P(IPD,5))**2)**3*SQRT(ABS(1.-REEM))*(1.+.5*REEM) IF(WTEE.LT.CJRN(0)) GO TO 201 PM(2,5)=EEM GO TO 260 ENDIF C 7 CALCULATION OF MAXIMUM WEIGHT WTMAX=1./WTCOR(ND) PMAX=PM(1,5)-PSUM+P(I+ND,5) PMIN=0. MLLW=1 MLUP=ND-1 DO 210 ML=MLLW,MLUP IL=MLUP-ML+1 PMAX=PMAX+P(I+IL,5) PMIN=PMIN+P(I+IL+1,5) 210 WTMAX=WTMAX*PAWT(PMAX,PMIN,P(I+IL,5)) C 8 M-GENERATOR GIVES WEIGHT, IF REJECTED TRY AGAIN 220 RND(1)=1. IL1LW=2 IL1UP=ND-1 DO 240 IL1=IL1LW,IL1UP RAN=CJRN(0) ML2LW=1 ML2UP=IL1-1 DO 230 ML2=ML2LW,ML2UP IL2=ML2UP-ML2+1 IF(RAN.LE.RND(IL2)) GOTO 240 230 RND(IL2+1)=RND(IL2) 240 RND(IL2+1)=RAN RND(ND)=0. WT=1. MLLW=1 MLUP=ND-1 DO 250 ML=MLLW,MLUP IL=MLUP-ML+1 PM(IL,5)=PM(IL+1,5)+P(I+IL,5)+(RND(IL)-RND(IL+1))*(P(IPD,5)-PSUM) 250 WT=WT*PAWT(PM(IL,5),PM(IL+1,5),P(I+IL,5)) IF(WT.LT.CJRN(0)*WTMAX) GOTO 220 C 9 TWO-PARTICLE DECAY IN CM 260 CONTINUE ILLW=1 ILUP=ND-1 DO 280 IL=ILLW,ILUP PA=PAWT(PM(IL,5),PM(IL+1,5),P(I+IL,5)) U(3)=2.*CJRN(0)-1. PHI=6.2832*CJRN(0) U(1)=SQRT(1.-U(3)**2)*COS(PHI) U(2)=SQRT(1.-U(3)**2)*SIN(PHI) DO 270 J=1,3 P(I+IL,J)=PA*U(J) 270 PM(IL+1,J)=-PA*U(J) P(I+IL,4)=SQRT(PA**2+P(I+IL,5)**2) 280 PM(IL+1,4)=SQRT(PA**2+PM(IL+1,5)**2) C 10 DECAY PRODUCTS LORENTZ TRANSFORMED TO LAB SYSTEM DO 290 J=1,4 290 P(I+ND,J)=PM(ND,J) MLLW=1 MLUP=ND-1 DO 320 ML=MLLW,MLUP IL=MLUP-ML+1 DO 300 J=1,3 300 BE(J)=PM(IL,J)/PM(IL,4) GA=PM(IL,4)/PM(IL,5) I1LW=I+IL I1UP=I+ND DO 320 I1=I1LW,I1UP BEP=BE(1)*P(I1,1)+BE(2)*P(I1,2)+BE(3)*P(I1,3) DO 310 J=1,3 310 P(I1,J)=P(I1,J)+GA*(GA/(1.+GA)*BEP+P(I1,4))*BE(J) 320 P(I1,4)=GA*(P(I1,4)+BEP) C 11 MATRIX ELEMENTS FOR SEMILEPTONIC AND OMEG AND PHI DECAYS C IN DECAY TABLE SET: NU 1ST IN C AND T DECAYS, CH. LEPTON 1ST IN B DECA IF(ND.EQ.3.AND.K(I+1,2).GE.7.AND.K(I+1,2).LE.18) THEN WT=FOUR(IPD,I+2)*FOUR(I+1,I+3) SW=P(I+1,5)**2+P(I+2,5)**2+2.*FOUR(I+1,I+2) SWX=P(IPD,5)**2 WT=WT*((SWX-PMAS(3)**2)**2+(PMAS(3)*WGAM)**2) + /((SW -PMAS(3)**2)**2+(PMAS(3)*WGAM)**2) IF(WT.LT.CJRN(0)*P(IPD,5)**4/WTCOR(1)) GOTO 220 ELSEIF(ND.EQ.3.AND.(K(IPD,2).EQ.92.OR.K(IPD,2).EQ.93)) THEN WT=(P(I+1,5)*P(I+2,5)*P(I+3,5))**2-(P(I+1,5)*FOUR(I+2,I+3))**2 + -(P(I+2,5)*FOUR(I+1,I+3))**2-(P(I+3,5)*FOUR(I+1,I+2))**2 + +2.*FOUR(I+1,I+2)*FOUR(I+1,I+3)*FOUR(I+2,I+3) IF(WT.LT.CJRN(0)*P(IPD,5)**6/WTCOR(2)) GOTO 220 ENDIF 330 CONTINUE IF(ND.EQ.0) RETURN DO 331 L=1,ND 331 K(I+L,1)=-IPD I=I+ND MNJTP=MAX(MNJTP,I) RETURN C C-- ABNORMAL EXIT 500 WRITE(ITLIS,501) MAXJTP,NEVENT 501 FORMAT(5(/),' NO. OF PARTICLES IN P( ,5) (/JET/) EXCEEDS',I10 1//1X,31HOVERFLOW OCCURRED IN SUB. DECAY 2//1X,11HEVENT NO. = ,I10 3//' INCREASE MAXJTP' 5) CALL OVERDM RETURN END