      PROGRAM SPMCGEN
C********************************************************************
C
C  Author: Stefan Ask, 2001
C
C  Toy Monte Carlo event generator for single photons from neutralino 
C  decay into a photon and a gravitino, hep-ph/9611437. 
C  e+e- --> XG --> GG + gamma  
C  X = Neutralino
C  G = Gravitino
C********************************************************************
      
      EXTERNAL GRNEXSEC
      EXTERNAL ONEFUNC
      EXTERNAL XSDIST

C Global variables -------------------------------
      REAL Sm,fmx,FMEL,FMER,xrfsusy
      REAL ECM,NEUM,SEL,SER,SQRTF
      REAL PI
      REAL TINT,ERR
      PARAMETER(PI=3.1415926536)
      PARAMETER(ECM=190.)
      PARAMETER(NEUM=100.)
      PARAMETER(SEL=300.)
      PARAMETER(SER=300.)
      PARAMETER(SQRTF=205.4)
      COMMON /SPVARS/ Sm,fmx,FMEL,FMER,xrfsusy
      COMMON /CINT/ TINT

C Random variables -------------------------------
C NTXRAN  = Neutralino theta distribution
C NPXRAN  = Neutralino phi distribution
C PTXRAN = Photon theta dist. in 
C          neutralino rest frame
C PPXRAN = Photon phi dist. in 
C          neutralino rest frame
C-------------------------------------------------
      INTEGER LEN
      PARAMETER(LEN=100000)
      REAL NTXRAN(LEN),NPXRAN(LEN)
      REAL PTXRAN(LEN),PPXRAN(LEN)
      REAL FSPACE(200)

C Photon variables in neutralino rest frame ------
      REAL EPN,APPN,PPN(3)
      REAL THEPN,PHIPN

C Neutralino variables in lab frame --------------
      REAL ENL,APNL,PNL(3)
      REAL THENL,PHINL

C Photon variables in lab frame ------
      REAL EPL,APPL,PPL(3)
      REAL THEPL,PHIPL

C Vectors for Lorenz transform
      REAL S(4),A(4),X(4)

C TEST variables --------------------------------- 
      LOGICAL TEST
      DATA TEST /.TRUE./
      REAL TTHETA,TCOSV,TDXSEC
      REAL PLTHE,CPLTHE,PLENE

C Hbook stuff. -----------------------------------
      INTEGER    NWPAW
      PARAMETER (NWPAW=1000000)
      REAL PAW
      COMMON / PAWC / PAW(NWPAW)
      
      CALL HLIMIT(NWPAW)
      CALL HDELET(0)
      CALL HBOOK1(1,'Photon Energy',60,0.,120.,0)
      CALL HBOOK1(2,'Photon theta angle, cos()',20,-1.,1.,0)
      CALL HBOOK1(3,'Isotropic photon dist. in restframe',10,0.,1.,0)
      CALL HBOOK1(4,'Neutralino theta angle, cos(), 
     +               generated',40,-1.,1.,0)
      CALL HBOOK1(5,'Neutralino theta angle, cos(), 
     +               from x-sec.',40,-1.,1.,0)
      CALL HBOOK1(6,'Neutralino theta angle, rad.',70,0.,7.,0)

C Init SUSY parameters. ---------------------------
      Sm = ECM**2
      fmx = NEUM 
      FMEL = SEL
      FMER = SER
      xrfsusy = SQRTF

C Integrate differential cross section for --------
C normalization, needs mathlib. -------------------
      CALL RADAPT(GRNEXSEC,-1.,1.,0,0.1,0.1,TINT,ERR)
      IF (TEST) THEN
         WRITE(*,*)TINT,ERR
      ENDIF

C Start generate random numbers, needs mathlib. ---
      CALL FUNLXP(XSDIST,FSPACE,-1.,1.)
      CALL FUNLUX(FSPACE,NTXRAN,LEN)
      CALL FUNLXP(ONEFUNC,FSPACE,0.,1.)
      CALL FUNLUX(FSPACE,NPXRAN,LEN)

      CALL FUNLXP(ONEFUNC,FSPACE,-1.,1.)
      CALL FUNLUX(FSPACE,PTXRAN,LEN)
      CALL FUNLXP(ONEFUNC,FSPACE,0.,1.)
      CALL FUNLUX(FSPACE,PPXRAN,LEN)

      DO I=1,LEN
         CALL HFILL(3,PTXRAN(I),0.,1.)
         CALL HFILL(4,NTXRAN(I),0.,1.)
      ENDDO

C Test for the distributions. ----------------------
      IF (TEST) THEN
         DO I=1,40 
            TCOSV = (REAL(I)-0.5)/20.-1.
            TDXSEC = GRNEXSEC(TCOSV)/TINT
            CALL HFILL(5,TCOSV,0.,TDXSEC)
         ENDDO
*         WRITE(*,*)TINT
*         WRITE(*,*)'N density at 0.:',GRNEXSEC(0.)/TINT
*         WRITE(*,*)'N density at -1.:',GRNEXSEC(-1.)/TINT
      ENDIF

C Start generate the events. -----------------------

      DO I=1,LEN
C Photon in NF. 
         EPN = NEUM/2.
         APPN = NEUM/2.
         THEPN = ACOS(PTXRAN(I))
         PHIPN = 2.*PI*PPXRAN(I)
         PPN(1) = APPN*SIN(THEPN)*SIN(PHIPN)
         PPN(2) = APPN*SIN(THEPN)*COS(PHIPN)
         PPN(3) = APPN*COS(THEPN)
C Neutralino in LF.
         ENL = ECM/2. + NEUM**2/(2.*ECM)
         APNL = ECM/2. - NEUM**2/(2.*ECM)
         THENL = ACOS(NTXRAN(I))
         PHINL = 2.*PI*NPXRAN(I)
         PNL(1) = APNL*SIN(THENL)*SIN(PHINL)
         PNL(2) = APNL*SIN(THENL)*COS(PHINL)
         PNL(3) = APNL*COS(THENL)
C Lorentz Boost it into LF (CM-frame)
         S(1) = PNL(1)
         S(2) = PNL(2)
         S(3) = PNL(3)
         S(4) = ENL
         A(1) = PPN(1)
         A(2) = PPN(2)
         A(3) = PPN(3)
         A(4) = EPN

         CALL SLOREN(NEUM,S,A,X)

C Make result histograms.
         IF (TEST) THEN
            PLTHE = ACOLL(X(1),X(2),X(3),0.,0.,1.,1)
            CPLTHE = COS(PLTHE)
            PLENE = X(4)
            CALL HFILL(1,PLENE,0.,1./100.)
            CALL HFILL(2,CPLTHE,0.,1./100.)
            CALL HFILL(6,THENL,0.,1./100.)
         ENDIF

      ENDDO

      IF (TEST) THEN
         CALL HRPUT(0,'spmcgen.hbook','N')
      ENDIF

      END

      REAL FUNCTION ACOLL(PX1,PY1,PZ1,PX2,PY2,PZ2,SIGN)
C***********************************************************************
C  Calculates the angle between two vectors, acollinearety 
C***********************************************************************

      INTEGER SIGN
      REAL PX1,PY1,PZ1
      REAL PX2,PY2,PZ2

      IF (SIGN.EQ.1) THEN
         ACOLL = ACOS((PX1*PX2+PY1*PY2+PZ1*PZ2)/
     +        (SQRT(PX1*PX1+PY1*PY1+PZ1*PZ1)
     +        *SQRT(PX2*PX2+PY2*PY2+PZ2*PZ2)))
      ELSE 
         ACOLL = ACOS(ABS(PX1*PX2+PY1*PY2+PZ1*PZ2)/
     +        (SQRT(PX1*PX1+PY1*PY1+PZ1*PZ1)
     +        *SQRT(PX2*PX2+PY2*PY2+PZ2*PZ2)))
      ENDIF

      RETURN
      END

      REAL FUNCTION GRNEXSEC(COSTHE)
C********************************************************************
C  Production cross section ds/dcostheta for neutralino-gravitino
C  From LNZ, hep-ph/9611437 
C  Taken from Susygen3/16
C********************************************************************
C
CIMP      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT NONE
*KEEP,SM.
      INTEGER i1
      REAL FMZ,GAMMAZ,SIN2W,ALPHA,E2,G2,PI,S,COSW,
     +     FLC,FRC
      REAL zr(4,4)
      REAL xrfsusy
      REAL fmel,fmer,fmx,sm,tm,costhe,um,sgs,sgt,sgu,dz2,redz,dselt1
     +     ,dsert1,dselu1,dseru1,dselt,dselu,dsert,dseru,e1,g1,fle,fre
     +     ,clr,xr,xl,fgg,fzz,fgz,ftt,fuu,fgt,fgu,fzt,fzu,fstu,conz
      COMMON /SPVARS/ Sm,fmx,FMEL,FMER,xrfsusy


C Set neutralino composition, photino.
        i1 = 1
*        zr(1,i1)=1.
*        zr(2,i1)=0.
C     ... a purely zino neutralino
*      zr(1,i1)=0.
*      zr(2,i1)=1.
C     ... a pure bino if minus is allowed ?
        zr(1,i1)=0.8663
        zr(2,i1)=-0.499

      grnexsec=0.

C Constants.
      FMZ=91.19
      GAMMAZ=2.497
      PI=3.1415926536
      ALPHA=1./128.
      E2=4.*PI*ALPHA
      SIN2W = .231243
      COSW = .876788
      G2=E2/SIN2W

C Selectron mass.
*      FMEL=150.
*      FMER=150.
C Neutralino mass
*      fmx=100.
C Susy breaking scale, sqrt(f) = sqrt(Mg*sqrt(3)*Mp)
C sqrt(f)=0.2054 => Mg = 1.E-5 eV
*      xrfsusy = 0.2054*1000
*      Sm=190.**2
      
      S=Sm
      if(sm.lt.fmx**2)return

      Tm=-0.5*(Sm-fmx**2)*(1.-costhe)
      Um=-0.5*(Sm-fmx**2)*(1.+costhe)

C
C phase space
C

      SGS=2.*sm*(sm-fmx**2)*(um**2+tm**2)
      SGT=1.*tm*(tm-fmx**2)*(tm**2      )
      SGU=1.*um*(um-fmx**2)*(um**2      )

C
C propagators
C

      DZ2=1./((S-FMZ**2)**2+FMZ**2*GAMMAZ**2)
      REDZ=(S-FMZ**2)*DZ2

      DSELT1=(TM-FMEL**2)
      DSERT1=(TM-FMER**2)
      DSELU1=(UM-FMEL**2)
      DSERU1=(UM-FMER**2)

      DSELT=1./DSELT1
      DSELU=1./DSELU1
      DSERT=1./DSERT1
      DSERU=1./DSERU1
C
C couplings
C
      e1=sqrt(e2)
      g1=sqrt(g2)/cosw
      FLE=-1./2.+SIN2W
      FRE=SIN2W
      CLR=FLE**2+FRE**2
      XR=+zr(1,i1)*e1-zr(2,i1)*g1*fre
      XL=-zr(1,i1)*e1+zr(2,i1)*g1*fle

C
C amplitudes
C


      fgg=+zr(1,i1)**2*e2*sgs/sm**2
C PPM bug, g2 should be g1**2 !! arggggghhhhh
C      fzz=+zr(2,i1)**2*g2*clr/2.d0*sgs*dz2
      fzz=+zr(2,i1)**2*g1**2*clr/2.*sgs*dz2
      fgz=-zr(1,i1)*zr(2,i1)*e1*g1*(fle+fre)*sgs*redz/s

      ftt=xr**2*sgt*dsert**2+xl**2*sgt*dselt**2
      fuu=xr**2*sgu*dseru**2+xl**2*sgu*dselu**2

      fgt=+zr(1,i1)*e1*(-tm)**3*2.*(xr*dsert-xl*dselt)
      fgu=+zr(1,i1)*e1*(-um)**3*2.*(xr*dseru-xl*dselu)

C PPM  bug ... the below terms miss out the propagator (1/(t-me**2))
C      fZt=-zr(2,i1)*g1*(-tm)**3*2.d0*redz*sm*(xr*fre-xl*fle)
C      fZu=-zr(2,i1)*g1*(-um)**3*2.d0*redz*sm*(xr*fre-xl*fle)
      fZt=-zr(2,i1)*g1*(-tm)**3*2.*redz*sm*(xr*fre*dsert-xl*fle*dselt)
      fZu=-zr(2,i1)*g1*(-um)**3*2.*redz*sm*(xr*fre*dseru-xl*fle*dselu)

      fstu=fgg+ftt+fuu+fzz+fgt+fgu+fzt+fzu+fgz


      CONZ=fstu*(s-fmx**2)/32./PI/S**2
C
C total
C
      grnexsec=conz*0.389E+9/xrfsusy**4/2.

      IF(grnexsec.LT.0.)grnexsec=0.

      RETURN
      END

      REAL FUNCTION ONEFUNC(X)
C********************************************************************
C     To generate an isotropic distribution with FUNLUX
C********************************************************************
      
      REAL X

      ONEFUNC = 1.
      
      RETURN
      END

      REAL FUNCTION XSDIST(X)
C********************************************************************
C     Differential cross section normalized to the total 
C     cross section.
C********************************************************************
      
      REAL X,TINT
      COMMON /CINT/ TINT

      XSDIST = GRNEXSEC(X)/TINT
      
      RETURN
      END

      SUBROUTINE SLOREN(NEUM,S,A,X)
C********************************************************************
C     Lorentz transform.      
C********************************************************************
      
      REAL NEUM,S(4),A(4),X(4)

      X(4) = (A(4)*S(4)+A(1)*S(1)+A(2)*S(2)+A(3)*S(3))/NEUM
      X(1) = A(1) + S(1)*(X(4)+A(4))/(S(4)+NEUM) 
      X(2) = A(2) + S(2)*(X(4)+A(4))/(S(4)+NEUM) 
      X(3) = A(3) + S(3)*(X(4)+A(4))/(S(4)+NEUM) 
      
      RETURN
      END
