SUBROUTINE SUSGEN( ierr ) C ---------------------------------------------------------------------- C Purpose: Generate an R-parity violating event. C C Notes: Binning in DIVON can cause events to be generated C in illegal kinematic regimes. These events are C trapped and re-generated without affecting the user. C The proton remnant is simulated by a diquark system C in the case of e-u events, and a proton plus quark in C the case of e-(sea quark) events. C C The photino is simulated in the common block C LUJETS by a (massive) neutrino. C C Errors: > 0 Fatal, end run. C < 0 Skip event. C C Author : JMB C Started 21-Feb-1991 C Last change 18-Feb-1992 C ---------------------------------------------------------------------- C IMPLICIT NONE C INCLUDE 'SUSANA.INC' INCLUDE 'PYTHIA57.INC' C INTEGER ierr, maxnum, i, seed, iway, c_type INTEGER part1, part2, part3 REAL dummy, randum, temp_xsec, temp_error REAL cx, cy, cz, mag, ulmass REAL p1(4), p2(4), p3(4), p4(4) REAL xy(2), rn32 C DOUBLE PRECISION xycheck(2), dfun, comp, xsec, gamma DOUBLE PRECISION xycheck(2), dfun, comp, gamma DOUBLE PRECISION pzp, pt, pzq, ephino DOUBLE PRECISION px, py, pz, pm, energy, theta, phi EXTERNAL ulmass C 1001 CONTINUE ierr = 0 C C - Increment event number event = event + 1 C C - Set function type. function_type = 1 C C -- If regenerate flag is set, DFUN must be repartitioned, as C a new bunch of decay energies has been generated. IF (regen) THEN maxnum = 100000 CALL PARTN( 2, low, high, 1.0, maxnum) CALL INTGRL( 2, 0, isustr(4), temp_xsec, temp_error ) regen = .FALSE. ENDIF C C -- Generate t and x of event. 1 CONTINUE CALL RANGEN(2,xy) xj = xy(1) tm = xy(2) xycheck(1) = xy(1) xycheck(2) = xy(2) C C -- Determine what type of event this is. dummy = dfun(2,xycheck) IF (dummy.LE.0.0) THEN C Rangen is being stupid - make it try again! GOTO 1 ENDIF dummy = RN32(randum)*(xseca+xsecb) C C Come here for electrons C IF (isustr(8). eq. 0) THEN IF (dummy .LT. xseca) THEN C This was an e- u type interaction. C Outgoing Quark type. IF ( MOD(ISUSTR(5),3).EQ.1 ) THEN outq_flav = 1 ELSE IF ( MOD(ISUSTR(5),3).EQ.2 ) THEN outq_flav = 3 ELSE IF ( MOD(ISUSTR(5),3).EQ.0 ) THEN outq_flav = 5 ENDIF ELSE C This was an e- dbar type interaction. IF ( ISUSTR(5).LE.3 ) THEN C Outgoing Quark type. outq_flav = -2 ELSE IF ( ISUSTR(5).LE.6 ) THEN C Outgoing Quark type. outq_flav = -4 ELSE IF ( ISUSTR(5).LE.9 ) THEN outq_flav = -6 ENDIF ENDIF ENDIF C C Come here for positrons C If (isustr(8). eq. 1) THEN IF (dummy .LT. xseca) THEN C This was an e+ d type interaction. C Outgoing Quark type. IF ( MOD(ISUSTR(5),3).EQ.1 ) THEN outq_flav = 2 ELSE IF ( MOD(ISUSTR(5),3).EQ.2 ) THEN outq_flav = 4 ELSE IF ( MOD(ISUSTR(5),3).EQ.0 ) THEN outq_flav = 6 ENDIF ELSE C This was an e+ ubar type interaction. IF ( ISUSTR(5).LE.3 ) THEN C Outgoing Quark type. outq_flav = -1 ELSE IF ( ISUSTR(5).LE.6 ) THEN C Outgoing Quark type. outq_flav = -3 ELSE IF ( ISUSTR(5).LE.9 ) THEN outq_flav = -5 ENDIF ENDIf ENDIF C Reset function type. Function_type = 0 C C Outgoing quark mass. outq_mass = ULMASS( outq_flav ) C C Store DIS event kinematics yj = -tm/(xj*epcms) q2j = -tm C C -- Add the remains of the proton to /LUJETS/. CALL SUSPEC( ierr ) IF (ierr.LT.0) THEN C Skip event GOTO 1001 ELSE IF (ierr.NE.0) THEN C End run GOTO 8900 ENDIF C C -- Calculate the momentum components of produced quark and photino, C distributing randomly in phi and adding to /LUJETS/. pzp = ( (msus(8)**2-tm)*xj*pp(4) - & 2.0*pl(4)*xj*(pl(4)*pp(4)-pl(3)*pp(3)) & - pl(4)*(tm-outq_mass**2) ) / & ( 2.0*xj*( pl(4)*pp(3)-pl(3)*pp(4) ) ) ephino = ( 2.0*pzp*pl(3) - tm + msus(8)**2 )/(2.0*pl(4)) IF ((ephino**2 - pzp**2 - msus(8)**2).LT.0.0) THEN IF (isustr(3).GT.0) THEN WRITE(isustr(1),*) 'SUSGEN4:Neg.Sqrt:Skipping an event' WRITE(isustr(1),*) ephino**2,'-',pzp**2,'-',msus(8)**2 WRITE(isustr(1),*) xj,tm ENDIF GOTO 1001 ENDIF pt = SQRT(ephino**2 - pzp**2 - msus(8)**2) pzq = pl(3) + xj*pp(3) - pzp C C Add the decay quark. k(n+1,1) = 1 k(n+1,2) = outq_flav phi = rn32(randum) * twopi p(n+1,1) = pt * cos(phi) p(n+1,2) = pt * sin(phi) p(n+1,3) = pzq p(n+1,4) = SQRT( pt**2 + pzq**2 + outq_mass**2 ) p(n+1,5) = outq_mass C Add the decay photino... (as a massive neutrino) K(n+2,1) = 1 k(n+2,2) = 12 p(n+2,1) = -p(n+1,1) p(n+2,2) = -p(n+1,2) p(n+2,3) = pzp p(n+2,4) = ephino p(n+2,5) = msus(8) n=n+2 C C Store photino 4-momentum n_photino = n Do 10, i = 1,4 pphino(i) = p(n,i) 10 CONTINUE C C -- Fragment the quarks and decay the particles (except the photino). CALL LUEXEC IF ((event.LT.2).AND.(isustr(3).GT.0)) CALL LULIST(1) C C Save the number of particles we have so far... First_inter = n C C -----------> Now propagate and decay the photino... C C - Calculate mean distance travelled before decay. gamma = pphino(4)/msus(8) decay_length = 2.998*gamma*lphino*1.0E+8 C C -- Calculate actual distance travelled before decay. CALL SUSDKL( IErr ) IF (ierr.GT.0) THEN C Reset error flag. IErr = 0 IF (isustr(3).GT.0) THEN WRITE(isustr(1),*) 'SUSDKL:Warning! Time out on decay' ENDIF ENDIF C C -- Calculate secondary vertex position. C Direction cosines mag = sqrt(pphino(1)**2 + pphino(2)**2 + pphino(3)**2) cx = pphino(1)/mag cy = pphino(2)/mag cz = pphino(3)/mag vphino(1) = distance * cx vphino(2) = distance * cy vphino(3) = distance * cz C C -- Now decide which particles the photino decayed into. CALL SUSBCH( 2, iway, ierr ) C IF ( (iway.EQ.1) .OR. (iway.EQ.3)) THEN c_type = -1 ELSE c_type = +1 ENDIF C o-------------------------o C | part1 = lepton | C | part2 = singlet quark | C | part3 = doublet quark | C o-------------------------o C IF (iway.LE.2) THEN C e part1 = 11*c_type IF ( MOD(ISUSTR(5),3).EQ.1 ) THEN C d part2 = -1*c_type ELSE IF ( MOD(ISUSTR(5),3).EQ.2 ) THEN C s part2 = -3*c_type ELSE IF ( MOD(ISUSTR(5),3).EQ.0 ) THEN C b part2 = -5*c_type ENDIF C IF ( ISUSTR(5).LE.3 ) THEN C ubar part3 = 2*c_type ELSE IF ( ISUSTR(5).LE.6 ) THEN C cbar part3 = 4*c_type ELSE IF ( ISUSTR(5).LE.9 ) THEN C tbar part3 = 6*c_type ENDIF C ELSE C C nu part1 = 12*c_type IF ( MOD(ISUSTR(5),3).EQ.1 ) THEN C d(bar) part2 = 1*c_type ELSE IF ( MOD(ISUSTR(5),3).EQ.2 ) THEN C s(bar) part2 = 3*c_type ELSE IF ( MOD(ISUSTR(5),3).EQ.0 ) THEN C t(bar) part2 = 5*c_type ENDIF C IF ( ISUSTR(5).LE.3 ) THEN C dbar part3 = -1*c_type ELSE IF ( ISUSTR(5).LE.6 ) THEN C cbar part3 = -3*c_type ELSE IF ( ISUSTR(5).LE.9 ) THEN C tbar part3 = -5*c_type ENDIF C ENDIF C C -- Calculate quark & positron momenta. m1 = ulmass(part1) m2 = ulmass(part2) m3 = ulmass(part3) C CALL SUSDKM( p1, p2, p3, pphino, iway, ierr) IF (IERR.LT.0) THEN IF (isustr(3).GT.0) THEN WRITE(isustr(1),*) 'Skipping an event' ENDIF GOTO 1001 ENDIF C C -- Flag photino as decayed... k(n_photino,1) = 11 C C -- Add quarks and positron to lujets. K(n+1,1) = 1 k(n+1,2) = part1 K(n+1,3) = n_photino p(n+1,1) = p1(1) p(n+1,2) = p1(2) p(n+1,3) = p1(3) p(n+1,4) = p1(4) p(n+1,5) = m1 C K(n+2,1) = 2 k(n+2,2) = part2 k(n+2,3) = n_photino p(n+2,1) = p2(1) p(n+2,2) = p2(2) p(n+2,3) = p2(3) p(n+2,4) = p2(4) p(n+2,5) = m2 C K(n+3,1) = 1 k(n+3,2) = part3 k(n+3,3) = n_photino p(n+3,1) = p3(1) p(n+3,2) = p3(2) p(n+3,3) = p3(3) p(n+3,4) = p3(4) p(n+3,5) = m3 C n = n+3 C -- Fragment these new quarks and decay particles. CALL LUEXEC C event_done = event_done + 1 C 8900 CONTINUE IF ((isustr(3).NE.0).OR.(event.GE.NEvent)) THEN CALL RN32OT(seed) WRITE(isustr(1),9001) seed ENDIF RETURN C 9001 FORMAT(8X,'SUSGEN: Random number seed = ',I10,' at the end' &, ' of this event') C END