DOUBLE PRECISION FUNCTION SUXSEC(N,xy) C --------------------------------------------------------------------------- C Purpose: Return the cross section for the R-parity violating process C a) e- q --> photino q', and C b) e- q'bar --> photino qbar, calculated at x, Q2 (input C arguments) and for HERA beam conditions. C C Notes: "a" diagrams interfere amongst themselves, as do the "b"s, C but "a"s don't interfere with "b"s. C C All incoming masses are set ~0. Outgoing particle masses are C accounted for (except spectator jet). C C Cross section = d(sigma)/(dx.dt), in inverse GeV. C C IErr = 0 --> All ok. C IErr = 1 --> Attempt to divide by zero. C C e- \ / photino C \ sq'(r) / C >---------->--------------------< C / \ C q / \ q' a1) C C e- --->--------------------------->---- photino C | C | se- C | C q --->--------------------------->---- q' a2) C C C e- --->--------------------------->---- q' C | C | sq C | C q --->--------------------------->---- photino a3) C C C C e- \ / photino C \ sqbar / C >---------->--------------------< C / \ C q'bar / \ qbar b1) C C e- --->--------------------------->---- photino C | C | se C | C q'bar --->--------------------------->---- qbar b2) C C C e- --->--------------------------->---- qbar C | C | sq' C | C q'bar --->--------------------------->---- photino b3) C C xq = Bjorken x variable. C t = -Q2 (cross channel) C C From parton cross section calculation by Herbi Dreiner, C Oxford Theoretical Physics C Author: JMB C Last change 20/2/92 C --------------------------------------------------------------------------- C IMPLICIT NONE C INCLUDE 'SUSANA.INC' C INTEGER i, ierr, N REAL density_1, density_2 REAL ulmass DOUBLE PRECISION xq, scale, xquark DOUBLE PRECISION Re, Rd, Rd2, Ru, Ru2 DOUBLE PRECISION vsum(4) DOUBLE PRECISION ephotino, pzp, pt DOUBLE PRECISION xy(2), terma, termb, termc DOUBLE PRECISION tplus, tminus DIMENSION xquark(-6:6) EXTERNAL LUDATA EXTERNAL SUSDAT EXTERNAL PYDATA C C - Load x and cos( theta ). xq = xy(1) tm = xy(2) C C - Check requested range. IF ( (xq.GT.high(1)) .OR. (xq.LT.low(1)) .OR. & (tm.GT.high(2)) .OR. (tm.LT.low(2)) ) THEN SUXSEC = 0.0 GOTO 8900 ENDIF IF (xq.EQ.1.0) xq = 0.99999 C C - Calculate parton 4-momentum. DO 1, i = 1,3 pq(i) = xq * pp(i) 1 CONTINUE pq(4) = sqrt( pq(1)**2 + pq(2)**2 + pq(3)**2 ) C C C - Calculate Mandelstam variables. C o-----------------------------------o C | tm = e-photino channel (=-Q^2) | C | sm = (e + q)^2 | C | um = (e - g) | C | Note that they are different for | C | the two processes (doublet or | C | singlet incoming) | C o-----------------------------------o DO 2, i=1,4 vsum(i) = pl(i) + pq(i) 2 CONTINUE sm = vsum(4)**2 - vsum(1)**2 - vsum(2)**2 - vsum(3)**2 um = mp2 + msin2 - sm - tm C C Error check. IF ( (abs(um-msus(ssq)**2).LT. Small) .OR. & (abs(tm-msus(7)**2).LT.Small)) THEN IErr = 1 WRITE(isustr(1),*) 'SUXSEC: Attempt to divide by zero!' GOTO 8900 ENDIF C C - Get the parton densities. scale = SQRT(sm) CALL PFTOPDG(xq, scale, xquark) C C For Electrons C IF(ISUSTR(8).EQ.0) THEN IF ( ISUSTR(5).LE.3 ) THEN density_2 = XQuark(2)/xq ELSE IF ( ISUSTR(5).LE.6 ) THEN density_2 = XQuark(4)/xq ELSE IF ( ISUSTR(5).LE.9 ) THEN density_2 = XQuark(6)/xq ENDIF IF ( MOD(ISUSTR(5),3).EQ.1 ) THEN density_1 = XQuark(-1)/xq ELSE IF ( MOD(ISUSTR(5),3).EQ.2 ) THEN density_1 = XQuark(-3)/xq ELSE IF ( MOD(ISUSTR(5),3).EQ.0 ) THEN density_1 = XQuark(-5)/xq ENDIF ENDIF C C For Positrons C IF(ISUSTR(8).EQ.1) THEN IF ( ISUSTR(5).LE.3 ) THEN density_2 = XQuark(1)/xq ELSE IF ( ISUSTR(5).LE.6 ) THEN density_2 = XQuark(3)/xq ELSE IF ( ISUSTR(5).LE.9 ) THEN density_2 = XQuark(5)/xq ENDIF IF ( MOD(ISUSTR(5),3).EQ.1 ) THEN density_1 = XQuark(-2)/xq ELSE IF ( MOD(ISUSTR(5),3).EQ.2 ) THEN density_1 = XQuark(-4)/xq ELSE IF ( MOD(ISUSTR(5),3).EQ.0 ) THEN density_1 = XQuark(-6)/xq ENDIF ENDIF C C -----> Calculate MATOT2 C - Check t+ & t-. terma = ( mp2 - msin2 )**2/(4.0*sm) termb = SQRT(sm/4.0) termc = (sm+mp2-msin2)**2/(4.0*sm) - mp2 IF ((termc.LT.0.0).OR.(wssq.EQ.0)) THEN XSECA = 0.0 MATOT2 = 0.0 ELSE tplus = terma - (termb + SQRT(termc))**2 tminus = terma - (termb - SQRT(termc))**2 IF ((tm.GT.tminus).OR.(tm.LT.tplus)) THEN XSECA = 0.0 MATOT2 = 0.0 ELSE Re = (tm-msus(7)**2) Ru = (um-msus(dsq2)**2) Rd2 = (sm-msus(ssq)**2)**2 + msus(ssq)**2 * wssq**2 MATOT2 = lambda**2*e2/2.0 * (ed**2*sm*(sm-mp2-msin2)/Rd2 & + (tm-msin2)*(tm-mp2)/Re**2 & + eu**2*(um-msin2)*(um-mp2)/Ru**2 & - 2.0*ed*ee*sm*tm*(sm-msus(ssq)**2)/(Rd2*Re) & - 2.0*ed*eu*sm*um*(sm-msus(ssq)**2)/(Rd2*Ru) & + 2.0*eu*ee*(tm*um-mp2*msin2)/(Re*Ru) ) C C - Calculate cross section. XSECA = 1.0/(16.0*pi*sm**2) * MATOT2 * density_2 * branch ENDIF ENDIF C C -- Crossing - and remember INCOMING QUARKS ONLY have zero mass. C Redefine for MBTOT2 um = mp2 + mdu2 - sm - tm C C - Check t+ & t-. terma = ( mp2 - mdu2 )**2/(4.0*sm) termc = (sm+mp2-mdu2)**2/(4.0*sm) - mp2 IF ((termc.LT.0.0).OR.(wdsq2.EQ.0.0)) THEN XSECB = 0.0 MBTOT2 = 0.0 ELSE tplus = terma - (termb + SQRT(termc))**2 tminus = terma - (termb - SQRT(termc))**2 IF ((tm.GT.tminus).OR.(tm.LT.tplus)) THEN C -- Do nothing XSECB = 0.0 MBTOT2 = 0.0 ELSE Re = (tm-msus(7)**2) Rd = (um-msus(ssq)**2) Ru2 = (sm-msus(dsq2)**2)**2 + msus(dsq2)**2 * wdsq2**2 MBTOT2 = lambda**2*e2/2.0 * (ed**2*(um-mdu2)*(um-mp2)/Rd**2 & + (tm-mdu2)*(tm-mp2)/Re**2 & + eu**2*sm*(sm-mp2-msin2)/Ru2 & - 2.0*ed*ee*(um*tm-mp2*mdu2)/(Rd*Re) & - 2.0*ed*eu*sm*um*(sm-msus(dsq2)**2)/(Ru2*Rd) & + 2.0*eu*ee*sm*tm*(sm-msus(dsq2)**2)/(Ru2*Re) ) C C - Calculate cross section. XSECB = 1.0/(16.0*pi*sm**2)*MBTOT2*density_1*branch ENDIF ENDIF C SUXSEC = XSECA + XSECB C C - Error checking. IF (MATOT2.lt.0.0) THEN WRITE (isustr(1),*) 'SUXSEC Error: MATOT2 =',MATOT2 WRITE (isustr(1),*) 'x, t',xq,tm ENDIF IF (MBTOT2.lt.0.0) THEN WRITE (isustr(1),*) 'SUXSEC Error: MBTOT2 =',MBTOT2 WRITE (isustr(1),*) 'x, t',xq,tm ENDIF C IF (isustr(3).GT.1) THEN write(6,*) 'Matrix Elements' write(6,*) matot2, mbtot2 write(6,*) 's,t,u' write(6,*) sm, tm, um write(6,*) 'Structure functions' write(6,*) density_1, density_2 ENDIF C 8900 CONTINUE C 9001 FORMAT( 2X,'SUXSEC:Error at x =',E15.4,' y =',E15.4,' MATOT2 = ' & , E15.5, ' Q2 =', E15.4) 9002 FORMAT( 2X,'SUXSEC:Error at x =',E15.4,' y =',E15.4,' MBTOT2 = ' & , E15.5, ' Q2 =', E15.4) 9003 FORMAT( 2X,'M1^2, M2^2, M3^2, M1M2, M1M3, M2M3',/, & E15.5,E15.5,E15.5,E15.5,E15.5,E15.5 ) C RETURN END