CDECK  ID>, UPINIT.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Chris Harris and Peter Richardson
*-- this version is modified so that it produces
*-- a eventually charged black hole remnant by
*-- Ben Koch 20/07/05
*-- Charge minimisation has been switched off.
*-- Changes and comments have been marked by the pattern **my*.

C----------------------------------------------------------------------
      SUBROUTINE UPINIT
C----------------------------------------------------------------------
C     Les Houches initialisation routine
C----------------------------------------------------------------------
      INCLUDE 'charybdis1001.inc'
C--Les Houches run common block
      INTEGER MAXPUP
      PARAMETER(MAXPUP=100)
      INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
      DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
      COMMON /HEPRUP/ IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
     &                IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),
     &                XMAXUP(MAXPUP),LPRUP(MAXPUP)
C--event common block
      INTEGER MAXNUP
      PARAMETER (MAXNUP=500)
      INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
      DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
      COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
     &              IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP),
     &              ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP),
     &              SPINUP(MAXNUP)
**my* Sabine Hossenfelders common block
      real*8 myM_f, myMrel
      common /constants/myM_f,myMrel
**endmy*
C--local variables
      INTEGER NSEARCH,I
      DOUBLE PRECISION WGTSUM,WGTSQ,SPNZRO,SPNHLF,SPNONE,PTOTAL
      DOUBLE PRECISION PPROBS(6,3)
      DATA PPROBS /1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0,
     &     0.773d0, 0.776d0, 0.761d0, 0.743d0, 0.726d0, 0.711d0,
     &     0.688d0, 0.827d0, 0.909d0, 0.959d0, 0.990d0, 1.011d0/
C--beam particles and energies (only proton=2212 and antiproton=-2212)
      IDBMUP(1) = 2212
      IDBMUP(2) = 2212
      EBMUP(1)  = 7000.0D0
      EBMUP(2)  = 7000.0D0
C--Set MPLNCK and define what is meant by it:
C--MSSDEF=1 means M_GT, MSSDEF=2 means M_DL and MSSDEF=3 means M_D
      MPLNCK=1000.0D0
**my**To asshure compability of the new energy spectrum of
**my**particles emitted during the Hawking phase with the
**my**rest of the charydis code MSSDEF has to be set to 1!
      MSSDEF=1
      myM_f=MPLNCK
      myMrel=2000.0d0
C--Set number of dimensions (number of EXTRA dimensions is n=D-4)
C--TOTDIM can be in the range 6-11
      TOTDIM=6
C--Use Giddings+Thomas momentum scale for calling pdfs?
C--(See page 12 of hep-ph/0106219)
      GTSCA=.FALSE.
C--Set mass window for black holes produced
      MINMSS=myMrel
      MAXMSS=13999.0D0
C--Set NBODY decay of BH remnant - NBODY can be in range 2-5.
      NBODY=2
C--Turn time variation of BH Hawking temperature in decay on or off
**my**the timevar=.FALSE. case doesn't really make sense for the
**my**simulation of remnant formation
      TIMVAR=.TRUE.
C--Set which decay products are to be used:
C--MSSDEC=1 gives no heavy particles
C--MSSDEC=2 gives t, W and Z as well
C--MSSDEC=3 gives t, W, Z and Higgs as well  
**my**MSSDEC is set dynamically to 2 during the program
**my**so this setting here is just dummy
      MSSDEC=2
C--Turn grey-body factors on/off
**my**GRYBDY factors have not been implemented
**my**so this variable should be false
      GRYBDY=.FALSE.
C--Turn kinematic cut-off of decay on (rather than M=MPLANCK cut-off)
      KINCUT=.FALSE.
C--print out option (0=no printout, 1 =errors only, 2=errors+info)
      IPRINT = 1
C--parameters for the initial maximum weight search
      NSEARCH=10000
C--pdf's for the beams (use the generator default)
C--MUST BE THE SAME FOR BOTH BEAMS
      PDFGUP(1) = -1
      PDFGUP(2) = -1
      PDFSUP(1) = -1
      PDFSUP(2) = -1
C--------------------------------------------------------------------
C     END OF USER VARIABLES DON'T TOUCH ANYTHING BELOW HERE
C--------------------------------------------------------------------
C--probabilities for particles of different spin (arbitrary normalisation)
      IF (GRYBDY) THEN
         SPNZRO = PPROBS(TOTDIM-5,1)
         SPNHLF = PPROBS(TOTDIM-5,2)
         SPNONE = PPROBS(TOTDIM-5,3)
      ELSE
         SPNZRO = 1d0
         SPNHLF = 0.750d0
         SPNONE = 1d0
      ENDIF
C--calculation of probability of emission for different particle types
C--only light SM particles
      IF (MSSDEC.EQ.1) THEN
         PTOTAL = (78.0d0*SPNHLF)+(18.0d0*SPNONE)
         PQUARK = (60.0d0*SPNHLF)/PTOTAL
         PLEPT  = (12.0d0*SPNHLF)/PTOTAL
         PNEUT  = (6.0d0*SPNHLF)/PTOTAL
         PGLUON = (16.0d0*SPNONE)/PTOTAL
         PGAMMA = (2.0d0*SPNONE)/PTOTAL
         PWBOSN = 0.0D0
         PZBOSN = 0.0D0
         PHIGGS = 0.0D0   
C--light SM partcles + top W/Z      
      ELSEIF (MSSDEC.EQ.2) THEN
         PTOTAL = (3.0d0*SPNZRO)+(90.0d0*SPNHLF)+(24.0d0*SPNONE)
         PQUARK = (72.0d0*SPNHLF)/PTOTAL
         PLEPT  = (12.0d0*SPNHLF)/PTOTAL
         PNEUT  = (6.0d0*SPNHLF)/PTOTAL
         PGLUON = (16.0d0*SPNONE)/PTOTAL
         PGAMMA = (2.0d0*SPNONE)/PTOTAL
         PZBOSN = (SPNZRO+(2.0d0*SPNONE))/PTOTAL
         PWBOSN = 2.0d0*PZBOSN
         PHIGGS = 0.0D0 
C--light SM particles +top W/Z and Higgs
      ELSE
         PTOTAL = (4.0d0*SPNZRO)+(90.0d0*SPNHLF)+(24.0d0*SPNONE)
         PQUARK = (72.0d0*SPNHLF)/PTOTAL
         PLEPT  = (12.0d0*SPNHLF)/PTOTAL
         PNEUT  = (6.0d0*SPNHLF)/PTOTAL
         PGLUON = (16.0d0*SPNONE)/PTOTAL
         PGAMMA = (2.0d0*SPNONE)/PTOTAL
         PZBOSN = (SPNZRO+(2.0d0*SPNONE))/PTOTAL
         PWBOSN = 2.0d0*PZBOSN
         PHIGGS = SPNZRO/PTOTAL
      ENDIF
C--what do do with the weights(here are generating weighted events)
      IDWTUP = 1
C--only one process
      NPRUP  = 1
C--communication code
      LPRUP(1) = 1
C--calculate the maximum weight
      WGTSUM = 0.0D0
      WGTSQ = 0.0D0
      XMAXUP(1) = 0.0D0
      DO I=1,NSEARCH
         CALL CHEVNT(.FALSE.)
         WGTSUM = WGTSUM+XWGTUP
         WGTSQ  = WGTSQ+XWGTUP**2
         XMAXUP(1) = MAX(XMAXUP(1),XWGTUP)
      ENDDO
      WGTSUM = WGTSUM/DBLE(NSEARCH)
      WGTSQ  = MAX(WGTSQ/DBLE(NSEARCH)-WGTSUM**2,0.0D0)
      WGTSQ  = SQRT(WGTSQ/ DBLE(NSEARCH))
C--cross section
      XSECUP(1) = wgtsum
C--error on the cross section
      XERRUP(1) = wgtsq
C--output initialisation information
C--header
      WRITE(*,10)
C--beam parameters
      WRITE(*,11) EBMUP(1),EBMUP(2)
C--basic parameters
      WRITE(*,12) MINMSS,MAXMSS,MPLNCK
     &,myM_f,myMrel
     &,MSSDEF,TOTDIM,TIMVAR,GRYBDY,
     &     KINCUT
C--choice of outgoing particles
      IF(MSSDEC.EQ.1) THEN
         WRITE(*,13)
      ELSEIF(MSSDEC.EQ.2) THEN
         WRITE(*,14)
      ELSEIF(MSSDEC.EQ.3) THEN
         WRITE(*,15)
      ENDIF
C--choice of scale
      IF(GTSCA) THEN
         WRITE(*,17)
      ELSE
         WRITE(*,18)
      ENDIF
C--information on the cross section
      WRITE(*,19) XSECUP(1),XERRUP(1),XMAXUP(1)
C--particle production probabilites
      WRITE(*,20) PQUARK,PLEPT,PNEUT,PGLUON,PGAMMA,PZBOSN,PWBOSN,PHIGGS
      RETURN
C--format statements for the output
 10   FORMAT(//10X,'       CHARYBDIS 1.000    ?? July 2003        ',//,
     &         10X,'Please reference: C.M. Harris, P. Richardson &',/,
     &         10X,'                  B.R. Webber                 ',/,
     &         10X,'                 hep-ph/0307???               '/)
 11   FORMAT(/10X,'INPUT CONDITIONS FOR THIS RUN'//
     &        10X,'BEAM1ENG. =',F10.2/
     &        10X,'BEAM2ENG. =',F10.2/)
 12   FORMAT(
     &     10X,'MINIMUM_BH_MASS    =',F11.2/
     &     10X,'MAXIMUM_BH_MASS    =',F11.2/
     &     10X,'PLANCK_MASS        =',F11.2/
     &     10X,'M_f        =',F11.2/
     &     10X,'M_Relic        =',F11.2/
     &     10X,'DEFN_OF_PLANCK_MASS =',I5/
     &     10X,'NO._OF_DIMENSIONS  =',I5/
     &     10X,'TIME_VARIATION     =',L5/
     &     10X,'GREY_BODY_EFFECTS  =',L5/
     &     10X,'KINEMATIC_CUT      =',L5)
 13   FORMAT(
     &     10X,'ONLY LIGHT SM PARTICLES PRODUCED')
 14   FORMAT(
     &     10X,'ALL_SM_PARTICLES_BUT_HIGGS_PRODUCED')
 15   FORMAT(
     &     10X,'ALL_SM_PARTICLES_PRODUCED')
      WRITE(*,16) NBODY
 16   FORMAT(
     &     10X,'PRODUCING ',I1,' PARTICLES IN REMNANT DECAY')
 17   FORMAT(
     &     10X,'USING_BLACK_RADIUS_AS_SCALE')
 18   FORMAT(
     &     10X,'USING_BLACK_MASS_AS_SCALE')
 19   FORMAT(
     &     10X,'CROSS_SECTION_(PB) =',G12.4/
     &     10X,'ERROR_IN_C-S _(PB) =',G12.4/
     &     10X,'MAXIMUM_WEIGHT     =',E12.4)
 20   FORMAT(/10X,'PARTICLE_PRODUCTION_PROBABILITIES'//
     &        10X,'QUARK                  =',F10.4/
     &        10X,'CHARGED_LEPTON         =',F10.4/
     &        10X,'NEUTRINO               =',F10.4/
     &        10X,'GLUON                  =',F10.4/
     &        10X,'PHOTON                 =',F10.4/
     &        10X,'Z_BOSON                =',F10.4/
     &        10X,'W_BOSON                =',F10.4/
     &        10X,'HIGGS BOSON            =',F10.4)
      END
CDECK  ID>, UPEVNT.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Chris Harris and Peter Richardson
C----------------------------------------------------------------------
      SUBROUTINE UPEVNT
C----------------------------------------------------------------------
C     Les Houches event routine
C----------------------------------------------------------------------
      IMPLICIT NONE
      CALL CHEVNT(.TRUE.)
      END
CDECK  ID>, CHDFIV.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Ian Knowles
C-----------------------------------------------------------------------
      SUBROUTINE CHDFIV(P0,P1,P2,P3,P4,P5,*)
C-----------------------------------------------------------------------
C     Generates 5-body decay 0->1+2+3+4+5 using pure phase space
C     (extracted from HERWIG)
C-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION CHRGEN,P0(5),P1(5),P2(5),P3(5),P4(5),P5(5),B,C,
     & AA,BB,CC,DD,EE,FF,TT,S1,RS1,GG,S2,RS2,HH,S3,PP,QQ,RR,SS,P1CM,
     & P2345(5),P2CM,P345(5),P3CM,P45(5),P4CM
      DOUBLE PRECISION TWO
      PARAMETER (TWO=2.D0)
      EXTERNAL CHRGEN
      B=P0(5)-P1(5)
      C=P2(5)+P3(5)+P4(5)+P5(5)
      IF (B.LT.C) RETURN 1
      AA=(P0(5)+P1(5))**2
      BB=B**2
      CC=C**2
      DD=(P3(5)+P4(5)+P5(5))**2
      EE=(P4(5)+P5(5))**2
      FF=(P4(5)-P5(5))**2
      TT=(B-C)*P0(5)**11/729
C Select squared masses S1, S2 and S3 of 2345, 345 and 45 subsystems
  10  S1=BB+CHRGEN(1)*(CC-BB)
      RS1=SQRT(S1)
      GG=(RS1-P2(5))**2
      S2=DD+CHRGEN(2)*(GG-DD)
      RS2=SQRT(S2)
      HH=(RS2-P3(5))**2
      S3=EE+CHRGEN(3)*(HH-EE)
      PP=(AA-S1)*(BB-S1)
      QQ=((RS1+P2(5))**2-S2)*(GG-S2)/S1
      RR=((RS2+P3(5))**2-S3)*(HH-S3)/S2
      SS=(S3-EE)*(S3-FF)/S3
      IF (PP*QQ*RR*SS*((GG-DD)*(HH-EE))**2.LT.TT*S1*S2*S3*CHRGEN(4)**2)
     & GOTO 10
C Do two body decays: 0-->1+2345, 2345-->2+345, 345-->3+45 and 45-->4+5
      P1CM=SQRT(PP/4)/P0(5)
      P2345(5)=RS1
      P2CM=SQRT(QQ/4)
      P345(5)=RS2
      P3CM=SQRT(RR/4)
      P45(5)=SQRT(S3)
      P4CM=SQRT(SS/4)
      CALL CHDTWO(P0   ,P1,P2345,P1CM,TWO,.TRUE.)
      CALL CHDTWO(P2345,P2,P345 ,P2CM,TWO,.TRUE.)
      CALL CHDTWO(P345 ,P3,P45  ,P3CM,TWO,.TRUE.)
      CALL CHDTWO(P45  ,P4,P5   ,P4CM,TWO,.TRUE.)
      END
CDECK  ID>, CHDFOR.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Ian Knowles
C-----------------------------------------------------------------------
      SUBROUTINE CHDFOR(P0,P1,P2,P3,P4,*)
C-----------------------------------------------------------------------
C     Generates 4-body decay 0->1+2+3+4 using pure phase space
C     (extracted from HERWIG)
C-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION CHRGEN,P0(5),P1(5),P2(5),P3(5),P4(5),B,C,AA,BB,
     & CC,DD,EE,TT,S1,RS1,FF,S2,PP,QQ,RR,P1CM,P234(5),P2CM,P34(5),P3CM
      DOUBLE PRECISION TWO
      PARAMETER (TWO=2.D0)
      EXTERNAL CHRGEN
      B=P0(5)-P1(5)
      C=P2(5)+P3(5)+P4(5)
      IF (B.LT.C) RETURN 1
      AA=(P0(5)+P1(5))**2
      BB=B**2
      CC=C**2
      DD=(P3(5)+P4(5))**2
      EE=(P3(5)-P4(5))**2
      TT=(B-C)*P0(5)**7/16
C Select squared masses S1 and S2 of 234 and 34 subsystems
  10  S1=BB+CHRGEN(1)*(CC-BB)
      RS1=SQRT(S1)
      FF=(RS1-P2(5))**2
      S2=DD+CHRGEN(2)*(FF-DD)
      PP=(AA-S1)*(BB-S1)
      QQ=((RS1+P2(5))**2-S2)*(FF-S2)/S1
      RR=(S2-DD)*(S2-EE)/S2
      IF (PP*QQ*RR*(FF-DD)**2.LT.TT*S1*S2*CHRGEN(3)**2) GOTO 10
C Do two body decays: 0-->1+234, 234-->2+34 and 34-->3+4
      P1CM=SQRT(PP/4)/P0(5)
      P234(5)=RS1
      P2CM=SQRT(QQ/4)
      P34(5)=SQRT(S2)
      P3CM=SQRT(RR/4)
      CALL CHDTWO(P0  ,P1,P234,P1CM,TWO,.TRUE.)
      CALL CHDTWO(P234,P2,P34 ,P2CM,TWO,.TRUE.)
      CALL CHDTWO(P34 ,P3,P4  ,P3CM,TWO,.TRUE.)
      END
CDECK  ID>, CHDTHR.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Bryan Webber
C-----------------------------------------------------------------------
      SUBROUTINE CHDTHR(P0,P1,P2,P3,*)
C-----------------------------------------------------------------------
C     GENERATES THREE-BODY DECAY 0->1+2+3 DISTRIBUTED
C     ACCORDING TO PHASE SPACE * WEIGHT
C     (extracted from HERWIG)
C-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION CHRGEN,CHRUNI,A,B,C,D,AA,BB,CC,DD,EE,FF,PP,QQ,WW,
     & RR,PCM1,PC23,P0(5),P1(5),P2(5),P3(5),P23(5),TWO,ONE
      EXTERNAL CHRGEN,CHRUNI
      PARAMETER (ONE=1.0D0,TWO=2.D0)
      A=P0(5)+P1(5)
      B=P0(5)-P1(5)
      C=P2(5)+P3(5)
      IF (B.LT.C) RETURN 1
      D=ABS(P2(5)-P3(5))
      AA=A*A
      BB=B*B
      CC=C*C
      DD=D*D
      EE=(B-C)*(A-D)
      A=0.5D0*(AA+BB)
      B=0.5D0*(CC+DD)
      C=4.0D0/(A-B)**2
C
C  CHOOSE MASS OF SUBSYSTEM 23 WITH PRESCRIBED DISTRIBUTION
C
   10 FF=CHRUNI(0,BB,CC)
      PP=(AA-FF)*(BB-FF)
      QQ=(CC-FF)*(DD-FF)
      WW=ONE
      RR=EE*FF*CHRGEN(0)
      IF (PP*QQ*WW.LT.RR*RR) GOTO 10
C
C  FF IS MASS SQUARED OF SUBSYSTEM 23.
C
C  DO 2-BODY DECAYS 0->1+23, 23->2+3
C
      P23(5)=SQRT(FF)
      PCM1=SQRT(PP)*0.5D0/P0(5)
      PC23=SQRT(QQ)*0.5D0/P23(5)
      CALL CHDTWO(P0,P1,P23,PCM1,TWO,.TRUE.)
      CALL CHDTWO(P23,P2,P3,PC23,TWO,.TRUE.)
      END
CDECK  ID>, CHDTWO.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Bryan Webber & Mike Seymour
C-----------------------------------------------------------------------
      SUBROUTINE CHDTWO(P0,P1,P2,PCM,COSTH,ZAXIS)
C-----------------------------------------------------------------------
C     GENERATES DECAY 0 -> 1+2
C
C     PCM IS CM MOMENTUM
C
C     COSTH = COS THETA IN P0 REST FRAME (>1 FOR ISOTROPIC)
C     IF ZAXIS=.TRUE., COS THETA IS MEASURED FROM THE ZAXIS
C     IF .FALSE., IT IS MEASURED FROM P0'S DIRECTION 
C     (extracted from HERWIG)
C-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION CHRUNI,ONE,ZERO,PCM,COSTH,C,S,P0(5),P1(5),P2(5),
     & PP(5),R(9)
      LOGICAL ZAXIS
      EXTERNAL CHRUNI
      PARAMETER (ZERO=0.D0, ONE=1.D0)
C--CHOOSE C.M. ANGLES
      C=COSTH
      IF (C.GT.ONE) C=CHRUNI(0,-ONE,ONE)
      S=SQRT(ONE-C*C)
      CALL CHRAZM(PCM*S,PP(1),PP(2))
C--PP IS MOMENTUM OF 2 IN C.M.
      PP(3)=-PCM*C
      PP(4)=SQRT(P2(5)**2+PCM**2)
      PP(5)=P2(5)
C--ROTATE IF NECESSARY
      IF (COSTH.LE.ONE.AND..NOT.ZAXIS) THEN
        CALL CHUROT(P0,ONE,ZERO,R)
        CALL CHUROB(R,PP,PP)
      ENDIF
C--BOOST FROM C.M. TO LAB FRAME
      CALL CHULOB(P0,PP,P2)
      CALL CHVDIF(4,P0,P2,P1)
      END
CDECK  ID>, CHEVNT.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Chris Harris
C----------------------------------------------------------------------
      SUBROUTINE CHEVNT(GENEVT)
C----------------------------------------------------------------------
C     Black Hole Production and Decay
C----------------------------------------------------------------------
      INCLUDE 'charybdis1001.inc'
C--Les Houches run common block
      INTEGER MAXPUP
      PARAMETER(MAXPUP=100)
      INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
      DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
      COMMON /HEPRUP/ IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
     &                IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),
     &                XMAXUP(MAXPUP),LPRUP(MAXPUP)
C--event common block
      INTEGER MAXNUP
      PARAMETER (MAXNUP=500)
      INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
      DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
      COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
     &              IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP),
     &              ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP),
     &              SPINUP(MAXNUP)
**my* Sabines common block
      real*8 myM_f, myMrel
      common /constants/myM_f,myMrel
C--Local Variables
      INTEGER myI
      DOUBLE PRECISION myKinCut
      DOUBLE PRECISION ECMS,ONE,TWO,PIFAC,THREE,ZERO,
     &     CHRGEN,RPOW,A0,A1,EMJ,FACT,QSQ,HCS,prob,
     &     FACTR,RCS,BHPOW,DISF1,DISF2,OMEGAD,DL2GT,
     &     RHORSQ,PLPOW,MPFACT,PA(5),PB(5),PC(5),PCM,CHUPCM,PCMF(5),
     &     XXMIN,XLMIN,EMSCA,XX(2),DISF(13,2),GEV2PB,CHFMAS
      INTEGER STAT,IBARYN,IQRK(MAXPUP),IGLU(MAXPUP),IQBR(MAXPUP),MTRY,
     &     NQRK,NGLU,NQBR,JQRK,JGLU,JQBR,J,ISTART,CHST,iline,FSTAT(5),
     &     IP,IQ,PT,I,BHCHRG,IDN(2),NTRY,NHTRY,SAVCHR,SAVBRN,ID(5),
     &     MHTRY,LTRY,LHTRY,CHFCHG,SPIN,FSPIN(5)
      LOGICAL GENEVT,CHRLOG,FIRST
      EXTERNAL CHFMAS,CHRGEN,CHUPCM,CHRLOG,CHFCHG
      PARAMETER(ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,ZERO=0.0D0,NHTRY=999)
      PARAMETER(MHTRY=20,LHTRY=200,GEV2PB=389379.D3)
      SAVE HCS,BHPOW,PLPOW,RPOW,A0,A1,FACT,MPFACT
      DATA FIRST/.TRUE./
      PIFAC = ACOS(-ONE)
      IF(FIRST) THEN
C--Convert MPLNCK to M_GT if above value is M_DL
         IF (MSSDEF.EQ.2) THEN
            DL2GT=(TWO**(TOTDIM-6)*PIFAC**(TOTDIM-5))**(1/(TOTDIM-TWO))
            INTMPL=DL2GT*MPLNCK
C--Convert MPLNCK to M_GT if above value is M_D
         ELSEIF (MSSDEF.EQ.3) THEN
            INTMPL=MPLNCK*TWO**(1/(TOTDIM-TWO))
         ELSE
            INTMPL=MPLNCK
         ENDIF
C--Calculate other parameters
         PLPOW=TWO/(TOTDIM-THREE)
         BHPOW=PLPOW-7.0D0
         RPOW=1/BHPOW
         A0=MINMSS**BHPOW
         A1=(MAXMSS**BHPOW-A0)
         IF (TOTDIM.EQ.11) THEN
            OMEGAD=(PIFAC**5)/12
         ELSEIF (TOTDIM.EQ.10) THEN
            OMEGAD=(32*PIFAC**4)/105
         ELSEIF (TOTDIM.EQ.9) THEN
            OMEGAD=(PIFAC**4)/3
         ELSEIF (TOTDIM.EQ.8) THEN
            OMEGAD=(16*PIFAC**3)/15
         ELSEIF (TOTDIM.EQ.7) THEN
            OMEGAD=PIFAC**3
         ELSEIF(TOTDIM.EQ.6) THEN
            OMEGAD=(8*PIFAC**2)/3
         ENDIF
         MPFACT=INTMPL**(TOTDIM-2)
         RHFACT=(4*(2*PIFAC)**(TOTDIM-4))/((TOTDIM-2)*OMEGAD*MPFACT)
         FIRST = .FALSE.
      ENDIF      
C--Select a mass for the produced black hole
      ECMS = 2.0D0*SQRT(EBMUP(1)*EBMUP(2))
      EMSCA=(A0+A1*CHRGEN(1))**RPOW
      BHMASS=EMSCA
      QSQ=EMSCA**2
      EMJ=EMSCA**(1-BHPOW)/BHPOW*A1
C--Calculate (r_h)**2 and initial T_H
      RHORSQ=(RHFACT*EMSCA)**PLPOW
C--Select initial momentum fractions
      XXMIN=QSQ/ECMS**2
      XLMIN=LOG(XXMIN)
      FACT=-GEV2PB*PIFAC*RHORSQ*EMJ*XLMIN*TWO/EMSCA
C--Change momentum scale for calling pdfs (see comment above)
      IF (GTSCA) EMSCA=1/SQRT(RHORSQ)
      SCALUP = EMSCA
      CALL CHPDF(XXMIN,XLMIN,EMSCA,XX,DISF)
      HCS=ZERO
      DISF1=ZERO
      DISF2=ZERO
      DO 30 IP=1,13
      DISF1=DISF1+DISF(IP,1)
 30   DISF2=DISF2+DISF(IP,2)
      FACTR=FACT*DISF1*DISF2
      HCS=FACTR
      XWGTUP = HCS
      IF(.NOT.GENEVT) RETURN
      RCS=HCS*CHRGEN(0)
      HCS = 0.0D0
      DO 20 IP=1,13
      DO 10 IQ=1,13
      FACTR=FACT*DISF(IP,1)*DISF(IQ,2)
      HCS=HCS+FACTR
      IF (HCS.GT.RCS) GOTO 99
 10   CONTINUE
 20   CONTINUE
C--Generate event
 99   IDN(1)=IP
      IDN(2)=IQ
C--put the incoming particles in the Les Houches common block
      NUP = 2
      DO I=1,2
         IF(IDN(I).GT.6.AND.IDN(I).LE.12) THEN
            IDN(I) = -MOD(IDN(I)-1,6)-1
         ELSEIF(IDN(I).EQ.13) THEN
            IDN(I) = 21
         ENDIF
         IDUP(I) = IDN(I)
         PUP(1,I) = 0.0D0
         PUP(2,I) = 0.0D0
         PUP(3,I) = XX(I)*EBMUP(I)
         PUP(4,I) = XX(I)*EBMUP(I)
         ISTUP(I) = -1
      ENDDO
      PUP(3,2) = -PUP(3,2)
      CALL CHVSUM(4,PUP(1,1),PUP(1,2),PCMF)
      CALL CHUMAS(PCMF)
C--Calculate BH charge (actually 3*charge)
      BHCHRG=CHFCHG(IDN(1))+CHFCHG(IDN(2))
C--Black Hole decay
**my***begin emitttion of Hawking radiation**************
***my**this KinCut forbids emission of particles
***with arbitray low energy, as it would happen in
***the end of the hawking phase.
      myKinCut=0.1d0
***my**MSSDEC=2 is without charge minimisation**
      MSSDEC=2
      MTRY = 0
      CHST = BHCHRG
 35   MTRY = MTRY+1
      BHCHRG = CHST
      NUP = 2
      CALL CHVEQU(5,PCMF,PA)
*       write(*,*) 'my 1 genBH-BHcharge_PA', BHCHRG, PA
 40   NUP = NUP+1
C--Check that BH mass > MPLNCK if not using KINCUT
      IF ((.NOT.KINCUT).AND.(PA(5).LT.(myMrel+myKinCut))
     &.AND.(BHCHRG/3d0.EQ.INT(BHCHRG/3d0))) then
*       Write(*,*) , 'my 2 End-of-Hawking-phase'
       Write(*,*) , 'my 3 BHCHRG-BHEnerg-BHMass-Mrel'
     & , BHCHRG, PA(4),PA(5),myMrel
       GOTO 55
      endif
C--Set BHMASS to correct value (necessary in time varying case)
      IF (TIMVAR) BHMASS=PA(5)
C--Select the decay product
      LTRY=0
 45   IF (MSSDEC.EQ.1) THEN
         CALL CHHBH1(PT,STAT,SPIN,BHCHRG)
      ELSEIF (MSSDEC.EQ.2) THEN
         CALL CHHBH2(PT,STAT,SPIN,BHCHRG)
      ELSE
         CALL CHHBH3(PT,STAT,SPIN,BHCHRG)
      ENDIF
C--Obtain energy of emitted parton
      PC(5)=CHFMAS(PT)
      CALL CHUBHS(PC(4),STAT,SPIN,PA(5),PC(5))
      If (PC(4).LT. myKinCut) then
       PC(4)=myKinCut
      ENDIF
**my* when emitted particle is a quark of light generation
**and its energy is smaller than 0.32 GeV so 
**take bare masses from PDG
      IF  ( (PC(4).LT.PC(5)) .AND. (Abs(PT).EQ.1)) THEN
       PC(5)=0.0029d0
      ELSEIF ((PC(4).LT.PC(5)).AND.(Abs(PT).EQ.2)) THEN
       PC(5)=0.0052d0
      ELSEIF ((PC(4).LT.PC(5)).AND.(Abs(PT).EQ.3)) THEN
       PC(5)=0.116d0
      ENDIF

C--Check energy > mass for selected parton
      IF (PC(4).LT.PC(5)) THEN
         BHCHRG=BHCHRG+CHFCHG(PT)
         IF(IPRINT.EQ.2) WRITE (*,*) 'Not enough energy'
         GOTO 45
      ENDIF
C--Check that emission is kinematically allowed
      IF (PC(4).GT.(PA(5)**2+PC(5)**2)/(TWO*PA(5))) THEN
         IF (KINCUT.OR.LTRY.GT.LHTRY) THEN 
            GOTO 50
         ELSE
            LTRY=LTRY+1
            BHCHRG=BHCHRG+CHFCHG(PT)
            GOTO 45
         ENDIF
      ENDIF
C--Use 2-body phase space decay 
      PB(5)=SQRT(PA(5)**2+PC(5)**2-TWO*PA(5)*PC(4))
C--Calculate centre yof mass momentum for decay
      PCM=SQRT(PC(4)**2-PC(5)**2)
C--PA is BH 5-momentum before, PB after decay
      CALL CHDTWO(PA,PB,PC,PCM,TWO,.TRUE.)
C--Add emitted particle to event record
      CALL CHVEQU(5,PC,PUP(1,NUP))
      IDUP(NUP)=PT
      ISTUP(NUP) = 1
      MOTHUP(1,NUP) = 1
      MOTHUP(2,NUP) = 2
C--Change BH 5-momentum back to PA
      CALL CHVEQU(5,PB,PA)
      GOTO 40
 50   BHCHRG=BHCHRG+CHFCHG(PT)
 55   NUP = NUP-1

***my*******the charge of the BH remnant is given
**my********to pythia by e+ (e-) which can be
**identified in the event analysis according to their
**unusual mass =0d0.
 57    If(BHCHRG.NE.0) then
       PC(1)=0d0
       PC(2)=0d0
       PC(3)=0d0
       PC(4)=0.001d0
       PC(5)=0.0d0
*CHFMAS(11)
       NUP=NUP+1
       PCM=SQRT(PC(4)**2-PC(5)**2)
**my**set New Momentum, Energy and Mass by hand
       DO I=1,4
         PB(I)=PA(I)-PC(I)
       ENDDO
       PB(5)=SQRT(PB(4)**2-PB(3)**2-PB(2)**2-PB(1)**2)
***
       CALL CHVEQU(5,PC,PUP(1,NUP))
       IDUP(NUP)=-11*SIGN(1,BHCHRG)
       ISTUP(NUP) = 1
       MOTHUP(1,NUP) = 1
       MOTHUP(2,NUP) = 2
       CALL CHVEQU(5,PB,PA)
       BHCHRG=BHCHRG-CHFCHG(IDUP(NUP))
      ENDIF
      If(BHCHRG.NE.0) goto 57
***********end emittion of Hawking phase*********************
*my******Put remnant to event record and
**name it a graviton!!
      DO I=1,3
       PC(I)=PA(I)
      ENDDO
      PC(5)=PA(5)-myKinCut
      PC(4)=SQRT(PC(1)**2+PC(2)**2+PC(3)**2+PC(5)**2)
      NUP=NUP+1
**my**set New Momentum, Energy and Mass by hand
      DO I=1,4
        PB(I)=PA(I)-PC(I)
      ENDDO
      PB(5)=SQRT(PB(4)**2-PB(3)**2-PB(2)**2-PB(1)**2)
***
      CALL CHVEQU(5,PC,PUP(1,NUP))
      IDUP(NUP)=39
      ISTUP(NUP) = 1
      MOTHUP(1,NUP) = 1
      MOTHUP(2,NUP) = 2
      CALL CHVEQU(5,PB,PA)
      BHCHRG=BHCHRG-0
******************************************

C--Calculate the baryon number violation
C--Find the baryon number of the inital state
      IBARYN = 0
      DO I=1,2
        IF(ABS(IDUP(I)).LE.6) IBARYN = IBARYN-SIGN(1,IDUP(I))
      ENDDO
C--and the final state
      DO I=3,NUP
        IF(ABS(IDUP(I)).LE.6) IBARYN = IBARYN+SIGN(1,IDUP(I))
      ENDDO
C--Try to balance this out
      NTRY = 0
      SAVCHR = BHCHRG
      SAVBRN = IBARYN
*my* if baryon Numer not conserved just do decay again.
      If(IBARYN.NE.0) then
**my*       write (*,*) 'my 01 Starting whole again'
       goto 35
      ENDIF
C--Now sort out the colour connections
C--First the initial state
      NQRK = 0
      NGLU = 0
      NQBR = 0
      DO I=1,2
        IF(IDUP(I).GT.0.AND.IDUP(I).LE.6) THEN
           NQBR = NQBR+1
           IQBR(NQBR) = I
        ELSEIF(IDUP(I).LT.0.AND.IDUP(I).GE.-6) THEN
           NQRK = NQRK+1
           IQRK(NQRK) = I
        ELSEIF(IDUP(I).EQ.21) THEN
           NGLU = NGLU+1
           IGLU(NGLU) = I
        ENDIF
      ENDDO
C--Then the final state
      DO I=3,NUP
        IF(IDUP(I).GT.0.AND.IDUP(I).LE.6) THEN
           NQRK = NQRK+1
           IQRK(NQRK) = I
        ELSEIF(IDUP(I).LT.0.AND.IDUP(I).GE.-6) THEN
           NQBR = NQBR+1
           IQBR(NQBR) = I
        ELSEIF(IDUP(I).EQ.21) THEN
           NGLU = NGLU+1
           IGLU(NGLU) = I
        ENDIF
      ENDDO
      
C--zero the colour connections
      DO I=1,NUP
         DO J=1,2
            ICOLUP(J,I) = 0
         ENDDO
      ENDDO
C--Now make the colour connections
      JGLU = 0
      JQBR = 0
      JQRK = 0
      IF(NQRK.GT.0) THEN 
         I = IQRK(1)
         JQRK = 1
         J = 1
      ELSEIF(NGLU.GT.0) THEN
         I = IGLU(1)
         JGLU = 1
         J = 2
      ELSE
         GOTO 200
      ENDIF
      ISTART = I
      ILINE  = 500
C--Find an antiquark or gluon to pair this with
 150   IF(J.EQ.1.OR.J.EQ.2) THEN
         IF(JQBR.EQ.NQBR.and.JGLU.EQ.NGLU) THEN
            GOTO 250
         ELSE
            PROB = DBLE(NQBR-JQBR)/DBLE(NQBR+NGLU-JQBR-JGLU)

C--Pair with an antiquark
            IF(CHRLOG(PROB)) THEN
               IF(JQBR.NE.NQBR) THEN
                 IF(NGLU.NE.JGLU.and.NQBR-JQBR.EQ.1) GOTO 150
                  JQBR = JQBR+1
                  IF(I.LE.2) THEN
                     ICOLUP(2,I         ) = ILINE
                  ELSE
                     ICOLUP(1,I         ) = ILINE
                  ENDIF
                  IF(IQBR(JQBR).LE.2) THEN
                     ICOLUP(1,IQBR(JQBR)) = ILINE
                  ELSE
                     ICOLUP(2,IQBR(JQBR)) = ILINE
                  ENDIF
                  ILINE = ILINE+1
                  I = IQBR(JQBR)
                  J = 3
                  GOTO 150
               ELSE
                  GOTO 250
               ENDIF
C--Pair with a gluon
            ELSE
               IF(JGLU.NE.NGLU) THEN
                  JGLU = JGLU+1
                  IF(I.LE.2) THEN
                     ICOLUP(2,I         ) = ILINE 
                  ELSE
                     ICOLUP(1,I         ) = ILINE 
                  ENDIF
                  IF(IGLU(JGLU).LE.2) THEN
                     ICOLUP(1,IGLU(JGLU)) = ILINE
                  ELSE
                     ICOLUP(2,IGLU(JGLU)) = ILINE
                  ENDIF
                  ILINE = ILINE+1
                  I = IGLU(JGLU)
                  J = 2
                  GOTO 150
               ELSE
                  GOTO 250
               ENDIF
            ENDIF
         ENDIF
C--Find a quark to pair this with
      ELSEIF(J.EQ.3) THEN
         IF(JQRK.NE.NQRK) THEN
            JQRK = JQRK+1
            IF(I.LE.2) THEN
               ICOLUP(2,I         ) = ILINE
            ELSE
               ICOLUP(1,I         ) = ILINE
            ENDIF
            IF(IQRK(JQRK).LE.2) THEN
               ICOLUP(1,IQRK(JQRK)) = ILINE
            ELSE
               ICOLUP(2,IQRK(JQRK)) = ILINE
            ENDIF
            ILINE = ILINE+1
            I = IQRK(JQRK)
            J = 1
            GOTO 150
         ELSE
            GOTO 250
         ENDIF
      ENDIF

 250  IF(I.LE.2) THEN
       
         ICOLUP(2,I)      = ILINE
      ELSE
         ICOLUP(1,I)      = ILINE
      ENDIF
      IF(ISTART.LE.2) THEN
         ICOLUP(1,ISTART) = ILINE
      ELSE
         ICOLUP(2,ISTART) = ILINE
      ENDIF
C--zero the flavours
 200  DO I=1,NUP
         IF(IDUP(I).GT.0.AND.IDUP(I).LE.6) THEN
            ICOLUP(2,I) = 0
         ELSEIF(IDUP(I).GE.-6.AND.IDUP(I).lt.0) THEN
            ICOLUP(1,I) = 0
         ELSEIF(IDUP(I).NE.21) THEN
            ICOLUP(1,I) = 0
            ICOLUP(2,I) = 0
         ENDIF
      ENDDO
      END
CDECK  ID>, CHFCHG.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Peter Richardson
C----------------------------------------------------------------------
      FUNCTION CHFCHG(ID)
C----------------------------------------------------------------------
C     Function to return the charge of a SM particle input is PDG code
C----------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER CHFCHG,ID
      IF(ABS(ID).LE.6) THEN
         IF(MOD(ID,2).EQ.0) THEN
            CHFCHG = 2
         ELSE
            CHFCHG = -1
         ENDIF
      ELSEIF(ABS(ID).GE.11.AND.ABS(ID).LE.16) THEN
         IF(MOD(ID,2).EQ.0) THEN
            CHFCHG = 0
         ELSE
            CHFCHG = -3
         ENDIF
      ELSEIF((ID.GE.21.AND.ID.LE.23).OR.ID.EQ.25) THEN
         CHFCHG = 0
      ELSEIF(ABS(ID).EQ.24) THEN
         CHFCHG = 3
      ELSE
         print *,'CHFCHG WARNING UNKNOWN PARTICLE',ID
         STOP
      ENDIF
      IF(ID.LT.0) CHFCHG = -CHFCHG
      END
CDECK  ID>, CHFMAS.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Peter Richardson
C----------------------------------------------------------------------
      FUNCTION CHFMAS(ID)
C----------------------------------------------------------------------
C     Function to return the mass of a SM particle input is PDG code
C     This is generator dependent
C----------------------------------------------------------------------
C--get the masses from the HERWIG common block
#if !defined(PYTHIA)
      INCLUDE 'HERWIG65.INC'
      INTEGER ID
      DOUBLE PRECISION CHFMAS
      IF(ABS(ID).LE.6) THEN
         CHFMAS = RMASS(ABS(ID))
      ELSEIF(ABS(ID).GE.11.AND.ABS(ID).LE.16) THEN
         CHFMAS = RMASS(110+ABS(ID))
      ELSEIF(ID.EQ.21) THEN
         CHFMAS = RMASS(13)
      ELSEIF(ID.EQ.22) THEN
         CHFMAS = 0.0D0
      ELSEIF(ID.EQ.23) THEN
         CHFMAS = RMASS(200)
      ELSEIF(ABS(ID).EQ.24) THEN
         CHFMAS = RMASS(198)
      ELSEIF(ID.EQ.25) THEN
         CHFMAS = RMASS(201)
      ELSE
         PRINT *,'CHFMAS WARNING UNKNOWN PARTICLE',ID
         STOP
      ENDIF
#endif
C--get the masses from the PYTHIA common block
#if defined(PYTHIA)
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)
      COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
      INTEGER PYCOMP,ID
      DOUBLE PRECISION CHFMAS
      EXTERNAL PYCOMP
      CHFMAS  = PMAS(PYCOMP(ID),1)
#endif
      END
CDECK  ID>, CHHBH1.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :   Chris Harris
C---------------------------------------------------------------------
      SUBROUTINE CHHBH1(PT,STAT,SPIN,BHCHRG)
C---------------------------------------------------------------------
C     Subroutine to select the type of the next particle
C---------------------------------------------------------------------
      INCLUDE 'charybdis1001.inc'
      INTEGER PT,STAT,SPIN,BHCHRG,CHFCHG,I
      DOUBLE PRECISION PTYPE,CHRGEN,EPS
      LOGICAL FIRST
      DATA FIRST/.TRUE./
      EXTERNAL CHRGEN,CHFCHG
      PARAMETER(EPS=1.0D-6)
C--initialisation on first call
      IF(FIRST) THEN
         PFERM (1) = PQUARK
         PFERM (2) = PFERM(1)+PLEPT
         PFERM (3) = PFERM(2)+PNEUT
         PBOSON(1) = PFERM(3)+PGLUON
         PBOSON(2) = PBOSON(1)+PGAMMA
C--check the sum
         IF(ABS(PBOSON(2)-1.0D0).GT.EPS) THEN
            IF(IPRINT.GE.1) THEN
               print *,'WARNING PROBABILITIES DO NOT SUM TO ONE'
               print *,'RENORMALISING'
            ENDIF
            DO I=1,3
               PFERM(I) = PFERM(I)/PBOSON(2)
            ENDDO
            PBOSON(1) = PBOSON(1)/PBOSON(2)
            PQUARK = PQUARK/PBOSON(2)
            PLEPT  = PLEPT /PBOSON(2)
            PNEUT  = PNEUT /PBOSON(2)
            PGLUON = PGLUON/PBOSON(2)
            PGAMMA = PGAMMA/PBOSON(2)
            PBOSON(2) = 1.0D0
         ENDIF
         FIRST = .FALSE.
      ENDIF
C--Assume decay product is a fermion
      STAT=-1
      SPIN=1
C--Choose decay product
      PTYPE=CHRGEN(2)
C--Change STAT to 1 and SPIN to 2 if we have a boson instead
      IF (PTYPE.GE.PFERM(3)) THEN 
         STAT=1
         SPIN=2
      ENDIF
C--Choose particle type
      IF (PTYPE.LT.PFERM(1)) THEN
C--Quarks
         PT=1+INT(5*PTYPE/PFERM(1))
C--Ensure BH charge remains close to 0, and antiparticles half time
         IF (((CHFCHG(PT)*BHCHRG).LT.0).OR.
     &       ((BHCHRG.EQ.0).AND.(CHRGEN(4).LT.0.5D0))) THEN 
            PT=-PT
         ENDIF
      ELSEIF (PTYPE.LT.PFERM(2)) THEN
C--Leptons
         PT=11+2*INT(3*(PTYPE-PFERM(1))/PLEPT)
C--Ensure BH charge remains close to 0, and antiparticles half time
         IF (((CHFCHG(PT)*BHCHRG).LT.0).OR.
     &       ((BHCHRG.EQ.0).AND.(CHRGEN(4).LT.0.5D0))) THEN 
            PT=-PT
         ENDIF
      ELSEIF (PTYPE.LT.PFERM(3)) THEN
C--Neutrinos
         PT=12+2*INT(3*(PTYPE-PFERM(2))/PNEUT)
         IF(CHRGEN(4).LT.0.5D0) PT = -PT
      ELSEIF (PTYPE.LT.PBOSON(1)) THEN
C--Gluons
         PT=21
      ELSE
C--Photons
         PT=22
      ENDIF
C--Calculate new charge of BH
      BHCHRG=BHCHRG-CHFCHG(PT)
      END
CDECK  ID>, CHHBH2.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :   Chris Harris
C---------------------------------------------------------------------
      SUBROUTINE CHHBH2(PT,STAT,SPIN,BHCHRG)
C---------------------------------------------------------------------
C     Subroutine to select the type of the next particle
C---------------------------------------------------------------------
      INCLUDE 'charybdis1001.inc'
      INTEGER PT,STAT,SPIN,BHCHRG,CHFCHG,I
      DOUBLE PRECISION PTYPE,CHRGEN,EPS
      LOGICAL FIRST
      DATA FIRST/.TRUE./
      EXTERNAL CHRGEN,CHFCHG
      PARAMETER(EPS=1.0D-6)
C--initialisation on first call
      IF(FIRST) THEN
         PFERM (1) = PQUARK
         PFERM (2) = PFERM(1)+PLEPT
         PFERM (3) = PFERM(2)+PNEUT
         PBOSON(1) = PFERM(3)+PGLUON
         PBOSON(2) = PBOSON(1)+PWBOSN
         PBOSON(3) = PBOSON(2)+PZBOSN
         PBOSON(4) = PBOSON(3)+PGAMMA
C--check the sum
         IF(ABS(PBOSON(4)-1.0D0).GT.EPS) THEN
            IF(IPRINT.GE.1) THEN
               print *,'WARNING PROBABILITIES DO NOT SUM TO ONE'
               print *,'RENORMALISING'
            ENDIF
            DO I=1,3
               PFERM(I) = PFERM(I)/PBOSON(4)
               PBOSON(I) = PBOSON(I)/PBOSON(4)
            ENDDO
            PQUARK = PQUARK/PBOSON(4)
            PLEPT  = PLEPT /PBOSON(4)
            PNEUT  = PNEUT /PBOSON(4)
            PGLUON = PGLUON/PBOSON(4)
            PGAMMA = PGAMMA/PBOSON(4)
            PZBOSN = PZBOSN/PBOSON(4)
            PWBOSN = PWBOSN/PBOSON(4)
            PBOSON(4) = 1.0D0
         ENDIF
         FIRST = .FALSE.
      ENDIF
C--Assume decay product is a fermion
      STAT=-1
      SPIN=1
C--Choose decay product
      PTYPE=CHRGEN(2)
C--Change STAT to 1 and SPIN to 2 if we have a boson instead except for one
C--degree of freedom of massive gauge bosons
      IF (PTYPE.GE.PFERM(3)) THEN 
         STAT=1
         SPIN=2
         IF ((PTYPE.GE.PBOSON(1)).AND.(PTYPE.LT.PBOSON(3))) THEN
            IF (3*CHRGEN(3).LT.1.0d0) SPIN=0
         ENDIF
      ENDIF
C--Choose particle type
      IF (PTYPE.LT.PFERM(1)) THEN
C--Quarks
         PT=1+INT(6*PTYPE/PFERM(1))
C--Ensure BH charge remains close to 0, and antiparticles half time

**my**Switch of Charge minimisation***
*((CHFCHG(PT)*BHCHRG).LT.0).OR.((BHCHRG.EQ.0).AND.  )
         IF (
     &       ((CHRGEN(4).LT.0.5D0))) THEN 
            PT=-PT
         ENDIF

      ELSEIF (PTYPE.LT.PFERM(2)) THEN
C--Leptons
         PT=11+2*INT(3*(PTYPE-PFERM(1))/PLEPT)
C--Ensure BH charge remains close to 0, and antiparticles half time

**my**Switch of Charge minimisation***
*((CHFCHG(PT)*BHCHRG).LT.0).OR.((BHCHRG.EQ.0).AND.
         IF (
     &       ((CHRGEN(4).LT.0.5D0))) THEN 
            PT=-PT
         ENDIF

      ELSEIF (PTYPE.LT.PFERM(3)) THEN
C--Neutrinos
         PT=12+2*INT(3*(PTYPE-PFERM(2))/PNEUT)
         IF(CHRGEN(4).LT.0.5D0) PT = -PT
      ELSEIF (PTYPE.LT.PBOSON(1)) THEN
C--Gluons
         PT=21
      ELSEIF (PTYPE.LT.PBOSON(2)) THEN
C--W+/W-
         PT=24
C--Ensure BH charge remains close to 0, and antiparticles half time

**my**Switch of Charge minimisation***
*((CHFCHG(PT)*BHCHRG).LT.0).OR.((BHCHRG.EQ.0).AND.
         IF (
     &       ((CHRGEN(4).LT.0.5D0))) THEN 
            PT=-PT
         ENDIF
      ELSEIF (PTYPE.LT.PBOSON(3)) THEN
C--Z0
         PT=23
      ELSE
C--Photons
         PT=22
      ENDIF
C--Calculate new charge of BH
      BHCHRG=BHCHRG-CHFCHG(PT)
      END
CDECK  ID>, CHHBH3.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :   Chris Harris
C---------------------------------------------------------------------
      SUBROUTINE CHHBH3(PT,STAT,SPIN,BHCHRG)
C---------------------------------------------------------------------
C     Subroutine to select the type of the next particle
C---------------------------------------------------------------------
      INCLUDE 'charybdis1001.inc'
      INTEGER PT,STAT,SPIN,BHCHRG,CHFCHG,I
      DOUBLE PRECISION PTYPE,CHRGEN,EPS
      LOGICAL FIRST
      DATA FIRST/.TRUE./
      EXTERNAL CHRGEN,CHFCHG
      PARAMETER(EPS=1.0D-6)
C--initialisation on first call
      IF(FIRST) THEN
         PFERM (1) = PQUARK
         PFERM (2) = PFERM(1)+PLEPT
         PFERM (3) = PFERM(2)+PNEUT
         PBOSON(1) = PFERM(3)+PGLUON
         PBOSON(2) = PBOSON(1)+PWBOSN
         PBOSON(3) = PBOSON(2)+PZBOSN
         PBOSON(4) = PBOSON(3)+PHIGGS
         PBOSON(5) = PBOSON(4)+PGAMMA
C--check the sum
         IF(ABS(PBOSON(5)-1.0D0).GT.EPS) THEN
            IF(IPRINT.GE.1) THEN
               print *,'WARNING PROBABILITIES DO NOT SUM TO ONE'
               print *,'RENORMALISING'
            ENDIF
            DO I=1,3
               PFERM(I) = PFERM(I)/PBOSON(5)
               PBOSON(I) = PBOSON(I)/PBOSON(5)
            ENDDO
            PBOSON(4) = PBOSON(4)/PBOSON(5)
            PQUARK = PQUARK/PBOSON(5)
            PLEPT  = PLEPT /PBOSON(5)
            PNEUT  = PNEUT /PBOSON(5)
            PGLUON = PGLUON/PBOSON(5)
            PGAMMA = PGAMMA/PBOSON(5)
            PZBOSN = PZBOSN/PBOSON(5)
            PWBOSN = PWBOSN/PBOSON(5)
            PHIGGS = PHIGGS/PBOSON(5)
            PBOSON(5) = 1.0D0
         ENDIF
         FIRST = .FALSE.
      ENDIF
C--Assume decay product is a fermion
      STAT=-1
      SPIN=1
C--Choose decay product
      PTYPE=CHRGEN(2)
C--Change STAT to 1 and SPIN to 2 if we have a boson instead except for one
C--degree of freedom of massive gauge bosons
      IF (PTYPE.GE.PFERM(3)) THEN
         STAT=1
         SPIN=2
         IF ((PTYPE.GE.PBOSON(1)).AND.(PTYPE.LT.PBOSON(3))) THEN
            IF (3*CHRGEN(3).LT.1.0d0) SPIN=0
         ENDIF
      ENDIF
C--Choose particle type
      IF (PTYPE.LT.PFERM(1)) THEN
C--Quarks
         PT=1+INT(6*PTYPE/PFERM(1))
C--Ensure BH charge remains close to 0, and antiparticles half time
         IF (((CHFCHG(PT)*BHCHRG).LT.0).OR.
     &       ((BHCHRG.EQ.0).AND.(CHRGEN(4).LT.0.5D0))) THEN 
            PT=-PT
         ENDIF
      ELSEIF (PTYPE.LT.PFERM(2)) THEN
C--Leptons
         PT=11+2*INT(3*(PTYPE-PFERM(1))/PLEPT)
C--Ensure BH charge remains close to 0, and antiparticles half time
         IF (((CHFCHG(PT)*BHCHRG).LT.0).OR.
     &       ((BHCHRG.EQ.0).AND.(CHRGEN(4).LT.0.5D0))) THEN 
            PT=-PT
         ENDIF
      ELSEIF (PTYPE.LT.PFERM(3)) THEN
C--Neutrinos
         PT=12+2*INT(3*(PTYPE-PFERM(2))/PNEUT)
         IF(CHRGEN(4).LT.0.5D0) PT = -PT
      ELSEIF (PTYPE.LT.PBOSON(1)) THEN
C--Gluons
         PT=21
      ELSEIF (PTYPE.LT.PBOSON(2)) THEN
C--W+/W-
         PT=24
C--Ensure BH charge remains close to 0, and antiparticles half time
         IF (((CHFCHG(PT)*BHCHRG).LT.0).OR.
     &       ((BHCHRG.EQ.0).AND.(CHRGEN(4).LT.0.5D0))) THEN 
            PT=-PT
         ENDIF
      ELSEIF (PTYPE.LT.PBOSON(3)) THEN
C--Z0
         PT=23
      ELSEIF (PTYPE.LT.PBOSON(4)) THEN
C--Higgs
         PT=25
         SPIN=0
      ELSE
C--Photons
         PT=22
      ENDIF
C--Calculate new charge of BH
      BHCHRG=BHCHRG-CHFCHG(PT)
      END
CDECK  ID>, CHPDF.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Peter Richardson
C----------------------------------------------------------------------
      SUBROUTINE CHPDF(XXMINB,XLMINB,EMSCAB,XXB,DISFB)
C----------------------------------------------------------------------
C     Subroutine to calculate the pdfs
C     generator dependent
C----------------------------------------------------------------------
#if !defined(PYTHIA)
      INCLUDE 'HERWIG65.INC'
      DOUBLE PRECISION EMSCAB,XXB(2),DISFB(13,2),XXMINB,XLMINB,
     &     CHRUNI,UPV,DNV,USEA,DSEA,STR,CHM,BTM,TOP,GLU
      INTEGER IDB(2),I,J
      CHARACTER *8 DUMMY
      DOUBLE PRECISION VALUE(20)
      CHARACTER*20 PARM(20)
      LOGICAL FIRST
      EXTERNAL CHRUNI
      DATA VALUE/20*0D0/,PARM/20*' '/
      DATA FIRST/.TRUE./
C--Les Houches run common block
      INTEGER MAXPUP
      PARAMETER(MAXPUP=100)
      INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
      DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
      COMMON /HEPRUP/ IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
     &                IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),
     &                XMAXUP(MAXPUP),LPRUP(MAXPUP)
C--use generator internal pdfs
      IF(PDFGUP(1).LT.0) THEN
C---SET UP INITIAL STATE
         CALL HWUIDT(1,IDBMUP(1),IDB(1),DUMMY)
         PART1=DUMMY
         CALL HWUIDT(1,IDBMUP(2),IDB(2),DUMMY)
         PART2=DUMMY
         IPART1 = IDB(1)
         IPART2 = IDB(2)
         NHEP=1
         ISTHEP(NHEP)=101
         PHEP(1,NHEP)=0.0D0
         PHEP(2,NHEP)=0.0D0
         PHEP(3,NHEP)=EBMUP(1)
         PHEP(4,NHEP)=EBMUP(1)
         PHEP(5,NHEP)=RMASS(IDB(1))
         JMOHEP(1,NHEP)=0
         JMOHEP(2,NHEP)=0
         JDAHEP(1,NHEP)=0
         JDAHEP(2,NHEP)=0
         IDHW(NHEP)=IPART1
         IDHEP(NHEP)=IDPDG(IDB(1))
         NHEP=NHEP+1
         ISTHEP(NHEP)=102
         PHEP(1,NHEP)=0.0D0
         PHEP(2,NHEP)=0.0D0
         PHEP(3,NHEP)=-EBMUP(2)
         PHEP(4,NHEP)=EBMUP(2)
         PHEP(5,NHEP)=RMASS(IDB(2))
         JMOHEP(1,NHEP)=0
         JMOHEP(2,NHEP)=0
         JDAHEP(1,NHEP)=0
         JDAHEP(2,NHEP)=0
         IDHW(NHEP)=IPART2
         IDHEP(NHEP)=IDPDG(IPART2)
C---NEXT ENTRY IS OVERALL CM FRAME
         NHEP=NHEP+1
         IDHW(NHEP)=14
         IDHEP(NHEP)=0
         ISTHEP(NHEP)=103
         JMOHEP(1,NHEP)=NHEP-2
         JMOHEP(2,NHEP)=NHEP-1
         JDAHEP(1,NHEP)=0
         JDAHEP(2,NHEP)=0
         CALL CHVSUM(4,PHEP(1,NHEP-1),PHEP(1,NHEP-2),PHEP(1,NHEP))
         CALL CHUMAS(PHEP(1,NHEP))
         XLMIN = XLMINB
         EMSCA = EMSCAB
         XXMIN = XXMINB
         CALL HWSGEN(.TRUE.)
         DO I=1,2
            XXB(I) = XX(I)
            DO J=1,13
               DISFB(J,I) = DISF(J,I)
            ENDDO
         ENDDO
      ELSE
C--initialisation
         IF(FIRST) THEN
            IF(PDFGUP(1).NE.PDFGUP(2).OR.PDFSUP(1).NE.PDFSUP(2)) THEN
               print *,'MUST HAVE SAME PDF FOR BOTH BEAMS'
               STOP
            ENDIF
            PARM(1)='NPTYPE'
            VALUE(1)=1
            PARM(2)='NGROUP'
            VALUE(2)=PDFGUP(1)
            PARM(3)='NSET'
            VALUE(3)=PDFSUP(1)
            CALL PDFSET(PARM,VALUE)
            FIRST = .FALSE.
         ENDIF
         XXB(1)=EXP(CHRUNI(0,0.0D0,XLMINB))
         XXB(2)=XXMINB/XXB(1)
         DO I=1,2
            CALL STRUCTM(XXB(I),EMSCAB,
     &           UPV,DNV,USEA,DSEA,STR,CHM,BTM,TOP,GLU)
            DISFB( 3,I) = STR
            DISFB( 4,I) = CHM
            DISFB( 5,I) = BTM
            DISFB( 6,I) = TOP
            DISFB( 9,I) = STR
            DISFB(10,I) = CHM
            DISFB(11,I) = BTM
            DISFB(12,I) = TOP
            DISFB(13,I) = GLU
            IF(IDBMUP(I).GT.0) THEN
               DISFB(1,I) = DNV+DSEA
               DISFB(2,I) = UPV+USEA
               DISFB(7,I) = DSEA
               DISFB(8,I) = USEA
            ELSE
               DISFB(1,I) = DSEA
               DISFB(2,I) = USEA
               DISFB(7,I) = DNV+DSEA
               DISFB(8,I) = UPV+USEA
            ENDIF
         ENDDO
      ENDIF
#endif
#if defined(PYTHIA)
      IMPLICIT NONE
C--Les Houches run common block
      INTEGER MAXPUP
      PARAMETER(MAXPUP=100)
      INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
      DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
      COMMON /HEPRUP/ IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
     &                IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),
     &                XMAXUP(MAXPUP),LPRUP(MAXPUP)
      DOUBLE PRECISION EMSCAB,XXB(2),DISFB(13,2),I,J,XXMINB,
     &     XLMINB,CHRUNI,XPQ(-25:25),Q2,UPV,DNV,USEA,DSEA,STR,CHM,BTM,
     &     TOP,GLU
      EXTERNAL CHRUNI
      DOUBLE PRECISION VALUE(20)
      CHARACTER*20 PARM(20)
      LOGICAL FIRST
      DATA VALUE/20*0D0/,PARM/20*' '/
      DATA FIRST/.TRUE./
      IF(PDFGUP(1).LT.0) THEN
      XXB(1)=EXP(CHRUNI(0,0.0D0,XLMINB))
      XXB(2)=XXMINB/XXB(1)
      Q2 = EMSCAB**2
      CALL PYPDFU(IDBMUP(1),XXB(1),Q2,XPQ)
      DO I=1,6
         DISFB(I,1) = XPQ(I)
         DISFB(I+6,1) = XPQ(I-6)
      ENDDO
      DISFB(13,1) = XPQ(21)
      CALL PYPDFU(IDBMUP(2),XXB(2),Q2,XPQ)
      DO I=1,6
         DISFB(I,2) = XPQ(I)
         DISFB(I+6,2) = XPQ(I-6)
      ENDDO
      DISFB(13,2) = XPQ(21)
      ELSE
C--initialisation
         IF(FIRST) THEN
            IF(PDFGUP(1).NE.PDFGUP(2).OR.PDFSUP(1).NE.PDFSUP(2)) THEN
               print *,'MUST HAVE SAME PDF FOR BOTH BEAMS'
               STOP
            ENDIF
            PARM(1)='NPTYPE'
            VALUE(1)=1
            PARM(2)='NGROUP'
            VALUE(2)=PDFGUP(1)
            PARM(3)='NSET'
            VALUE(3)=PDFSUP(1)
            CALL PDFSET(PARM,VALUE)
            FIRST = .FALSE.
         ENDIF
         XXB(1)=EXP(CHRUNI(0,0.0D0,XLMINB))
         XXB(2)=XXMINB/XXB(1)
         DO I=1,2
            CALL STRUCTM(XXB(I),EMSCAB,
     &           UPV,DNV,USEA,DSEA,STR,CHM,BTM,TOP,GLU)
            DISFB( 3,I) = STR
            DISFB( 4,I) = CHM
            DISFB( 5,I) = BTM
            DISFB( 6,I) = TOP
            DISFB( 9,I) = STR
            DISFB(10,I) = CHM
            DISFB(11,I) = BTM
            DISFB(12,I) = TOP
            DISFB(13,I) = GLU
            IF(IDBMUP(I).GT.0) THEN
               DISFB(1,I) = DNV+DSEA
               DISFB(2,I) = UPV+USEA
               DISFB(7,I) = DSEA
               DISFB(8,I) = USEA
            ELSE
               DISFB(1,I) = DSEA
               DISFB(2,I) = USEA
               DISFB(7,I) = DNV+DSEA
               DISFB(8,I) = UPV+USEA
            ENDIF
         ENDDO
      ENDIF
#endif
      END
CDECK  ID>, CHRAZM.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Bryan Webber
C-----------------------------------------------------------------------
      SUBROUTINE CHRAZM(PT,PX,PY)
C-----------------------------------------------------------------------
C     RANDOMLY ROTATED 2-VECTOR (PX,PY) OF LENGTH PT
C     (taken from HERWIG)
C-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION CHRGEN,PT,PX,PY,C,S,CS,QT,ONE,ZERO
      PARAMETER(ONE=1.0D0, ZERO=0.0D0)
      EXTERNAL CHRGEN
   10 C=2.0D0*CHRGEN(1)-1.0D0
      S=2.0D0*CHRGEN(2)-1.0D0
      CS=C*C+S*S
      IF (CS.GT.ONE .OR. CS.EQ.ZERO) GOTO 10
      QT=PT/CS
      PX=(C*C-S*S)*QT
      PY=2.0D0*C*S*QT
      END
CDECK  ID>, CHRGEN.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Peter Richardson
C-----------------------------------------------------------------------
      FUNCTION CHRGEN(I)
C-----------------------------------------------------------------------
C     random number generator (just calls HERWIG or PYTHIA one)
C-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION CHRGEN
      INTEGER I
C--use HERWIG random number generator
#if !defined(PYTHIA)
      DOUBLE PRECISION HWRGEN
      EXTERNAL HWRGEN
      CHRGEN = HWRGEN(I)
#endif
C--use the PYTHIA one
#if defined(PYTHIA)
      DOUBLE PRECISION PYR
      EXTERNAL PYR
      CHRGEN =  PYR(I)
#endif
      END
CDECK  ID>, CHRLOG.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Bryan Webber
C-----------------------------------------------------------------------
      FUNCTION CHRLOG(A)
C-----------------------------------------------------------------------
C     Returns .TRUE. with probability A
C     (taken from HERWIG)
C-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION CHRGEN,A,R
      LOGICAL CHRLOG
      EXTERNAL CHRGEN
      CHRLOG=.TRUE.
      R=CHRGEN(0)
      IF(R.GT.A) CHRLOG=.FALSE.
      END
CDECK  ID>, CHRUNI.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Bryan Webber
C-----------------------------------------------------------------------
      FUNCTION CHRUNI(I,A,B)
C-----------------------------------------------------------------------
C     Uniform random random number in range [A,B]
C     (taken from HERWIG)
C-----------------------------------------------------------------------
      DOUBLE PRECISION CHRUNI,CHRGEN,A,B,RN
      INTEGER I
      EXTERNAL CHRGEN
      RN=CHRGEN(I)
      CHRUNI=A+RN*(B-A)
      END
CDECK  ID>, CHUBHS.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :   Chris Harris
C---------------------------------------------------------------------
      SUBROUTINE CHUBHS(ENERGY,STAT,SPIN,myBHmass,myPC5)
C---------------------------------------------------------------------
C     Calculates energy spectrum for a black hole of given radius
C---------------------------------------------------------------------
      INCLUDE 'charybdis1001.inc'
      DOUBLE PRECISION ENERGY,RADHOR,HWKTMP,ETRAT
      DOUBLE PRECISION CHRGEN,CHRUNI,CHUTAB
      DOUBLE PRECISION mySPCVAL, myBHmass, mySPCMAX, myMXARRY(7)
      DOUBLE PRECISION Specf, myMrel, myM_f,myPC5,ENERGYMAX
      INTEGER myNTRY
      INTEGER NARRAY(6),STAT,SPIN,POINTS,I
      EXTERNAL CHRGEN,CHRUNI,CHUTAB, Specf
      common /constants/ myM_f, myMrel
      DOUBLE PRECISION ONE,TWO,THREE,ZERO,PIFAC,HALF
      PARAMETER(ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,ZERO=0.0D0,HALF=0.5D0)
      DATA NARRAY /151,121,131,141,119,151/
*my** these maximum values have been found by scanning over 
*M_BH and Energy for myM_f=2TeV and M_rel=1TeV
*for lower myM_f (or higer M_rel) lower myMXARRY values
*would be sufficent.

      myMXARRY(1)=0d0
      myMXARRY(2)=363149.979d0
      myMXARRY(3)=766883.929d0
      myMXARRY(4)=1195483.87d0
      myMXARRY(5)=1490075.17d0
      myMXARRY(6)=1551143.44d0
      myMXARRY(7)=1385850.23d0
      myNTRY=0
      PIFAC = ACOS(-ONE)

C     Spectrum
      RADHOR = (RHFACT*BHMASS)**(ONE/(TOTDIM-THREE))

      HWKTMP = (TOTDIM-THREE)/(4.0D0*PIFAC*RADHOR)
      IF(GRYBDY) THEN
      ELSE
         IF (STAT.EQ.-1) THEN
         ELSE
         ENDIF
         mySPCMAX=myMXARRY(TOTDIM-4)
      ENDIF
      POINTS=NARRAY(TOTDIM-5)
C     MC for energy
      ENERGYMAX=((myBHmass**2+myPC5**2-myMrel**2)/(2*myBHmass))
      If (ENERGYMAX.lt.0d0) then
       myNTRY=20000
       GOTO 21
      ENDIF
 10    myNTRY=myNTRY+1
**my**Energy Spectrum due to Sabine Hossenfelder***
**my neglect Mass of emitted particle
      ENERGY=((myBHmass**2+myPC5**2-myMrel**2)/(2*myBHmass))*CHRGEN(1)
      IF (GRYBDY) THEN
      ELSE
         mySPCVAL=Specf(TOTDIM-4,myBHmass,Energy,1d0,RADHOR)
      ENDIF

      IF (Abs((myBHmass-myMrel)/myMrel) .LT. 0.1d0) then
       IF (
     & (mySPCVAL.LT.Specf(TOTDIM-4,myBHmass,ENERGYMAX,1d0,RADHOR)
     &*CHRGEN(1))
     &  .AND.(myNTRY.LT.10000)) GOTO 10
       ELSE
        IF ((mySPCVAL.LT.mySPCMAX*CHRGEN(1))
     &  .AND.(myNTRY.LT.10000)) GOTO 10
      ENDIF
 21   If (myNTRY.GT.9999) then
**my**for numerics Reasons take energy to be very small
** if it was not possible in 9999 tries to find the a right one
       ENERGY=0.001d0
       If (myNTRY.ne.20000) then
*        write (*,*) , 'my 0 CHUBHS hat nicht genug gewruerfelt'
       endif
      ENDIF
      END
CDECK  ID>, CHULB4.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Adapted by Bryan Webber
C-----------------------------------------------------------------------
      SUBROUTINE CHULB4(PS,PI,PF)
C-----------------------------------------------------------------------
C     TRANSFORMS PI (GIVEN IN REST FRAME OF PS) INTO PF (IN LAB)
C     N.B. P(1,2,3,4) = (PX,PY,PZ,E); PS(5)=M
C     (taken from HERWIG)
C-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION PF4,FN,PS(5),PI(4),PF(4)
      IF (PS(4).EQ.PS(5)) THEN
        PF(1)= PI(1)
        PF(2)= PI(2)
        PF(3)= PI(3)
        PF(4)= PI(4)
      ELSE
        PF4  = (PI(1)*PS(1)+PI(2)*PS(2)
     &         +PI(3)*PS(3)+PI(4)*PS(4))/PS(5)
        FN   = (PF4+PI(4)) / (PS(4)+PS(5))
        PF(1)= PI(1) + FN*PS(1)
        PF(2)= PI(2) + FN*PS(2)
        PF(3)= PI(3) + FN*PS(3)
        PF(4)= PF4
      END IF
      END
CDECK  ID>, CHULOB.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Adapted by Bryan Webber
C-----------------------------------------------------------------------
      SUBROUTINE CHULOB(PS,PI,PF)
C-----------------------------------------------------------------------
C     TRANSFORMS PI (GIVEN IN REST FRAME OF PS) INTO PF (IN LAB)
C     N.B. P(1,2,3,4,5) = (PX,PY,PZ,E,M)
C     (taken from HERWIG)
C-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION PS(5),PI(5),PF(5)
      CALL CHULB4(PS,PI,PF)
      PF(5)= PI(5)
      END
CDECK  ID>, CHUMAS.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Bryan Webber
C-----------------------------------------------------------------------
      SUBROUTINE CHUMAS(P)
C-----------------------------------------------------------------------
C     PUTS INVARIANT MASS IN 5TH COMPONENT OF VECTOR
C     (NEGATIVE SIGN IF SPACELIKE) (extracted from HERWIG)
C-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION CHUSQR,P(5)
      EXTERNAL CHUSQR
      P(5)=CHUSQR((P(4)+P(3))*(P(4)-P(3))-P(1)**2-P(2)**2)
      END
CDECK  ID>, CHUPCM.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Bryan Webber
C-----------------------------------------------------------------------
      FUNCTION CHUPCM(EM0,EM1,EM2)
C-----------------------------------------------------------------------
C     C.M. MOMENTUM FOR DECAY MASSES EM0 -> EM1 + EM2
C     SET TO -1 BELOW THRESHOLD (extracted from HERWIG)
C-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION CHUPCM,EM0,EM1,EM2,EMS,EMD
      EMS=ABS(EM1+EM2)
      EMD=ABS(EM1-EM2)
      IF (EM0.LT.EMS.OR.EM0.LT.EMD) THEN
        CHUPCM=-1.0D0
      ELSEIF (EM0.EQ.EMS.OR.EM0.EQ.EMD) THEN
        CHUPCM=0.0D0
      ELSE
        CHUPCM=SQRT((EM0+EMD)*(EM0-EMD)*
     &              (EM0+EMS)*(EM0-EMS))*0.5D0/EM0
      ENDIF
      END
CDECK  ID>, CHUROB.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Bryan Webber
C-----------------------------------------------------------------------
      SUBROUTINE CHUROB(R,P,Q)
C-----------------------------------------------------------------------
C     ROTATES VECTORS BY INVERSE OF ROTATION MATRIX R 
C     (taken from HERWIG)
C-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION S1,S2,S3,R(3,3),P(3),Q(3)
      S1=P(1)*R(1,1)+P(2)*R(2,1)+P(3)*R(3,1)
      S2=P(1)*R(1,2)+P(2)*R(2,2)+P(3)*R(3,2)
      S3=P(1)*R(1,3)+P(2)*R(2,3)+P(3)*R(3,3)
      Q(1)=S1
      Q(2)=S2
      Q(3)=S3
      END
CDECK  ID>, CHUROT.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Bryan Webber
C-----------------------------------------------------------------------
      SUBROUTINE CHUROT(P,CP,SP,R)
C-----------------------------------------------------------------------
C     R IS ROTATION MATRIX TO GET FROM VECTOR P TO Z AXIS, FOLLOWED BY
C     A ROTATION BY PSI ABOUT Z AXIS, WHERE CP = COS-PSI, SP = SIN-PSI
C     (extracted from HERWIG)
C-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION WN,CP,SP,PTCUT,PP,PT,CT,ST,CF,SF,P(3),R(3,3)
      DATA WN,PTCUT/1.D0,1.D-20/
      PT=P(1)**2+P(2)**2
      PP=P(3)**2+PT
      IF (PT.LE.PP*PTCUT) THEN
         CT=SIGN(WN,P(3))
         ST=0.0D0
         CF=1.0D0
         SF=0.0D0
      ELSE
         PP=SQRT(PP)
         PT=SQRT(PT)
         CT=P(3)/PP
         ST=PT/PP
         CF=P(1)/PT
         SF=P(2)/PT
      END IF
      R(1,1)= CP*CF*CT+SP*SF
      R(1,2)= CP*SF*CT-SP*CF
      R(1,3)=-CP*ST
      R(2,1)=-CP*SF+SP*CF*CT
      R(2,2)= CP*CF+SP*SF*CT
      R(2,3)=-SP*ST
      R(3,1)= CF*ST
      R(3,2)= SF*ST
      R(3,3)= CT
      END
CDECK  ID>, CHUSQR.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Bryan Webber
C-----------------------------------------------------------------------
      FUNCTION CHUSQR(X)
C-----------------------------------------------------------------------
C     SQUARE ROOT WITH SIGN RETENTION (extracted from HERWIG)
C-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION CHUSQR,X
      CHUSQR=SIGN(SQRT(ABS(X)),X)
      END
CDECK  ID>, CHUTAB.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Adapted by Bryan Webber
C-----------------------------------------------------------------------
      FUNCTION CHUTAB(F,A,NN,X,MM)
C-----------------------------------------------------------------------
C     MODIFIED CERN INTERPOLATION ROUTINE DIVDIF
C-----------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER NN,MM,MMAX,N,M,MPLUS,IX,IY,MID,NPTS,IP,I,J,L,ISUB
      DOUBLE PRECISION CHUTAB,SUM,X,F(NN),A(NN),T(20),D(20)
      LOGICAL EXTRA
      DATA MMAX/10/
      N=NN
      M=MIN(MM,MMAX,N-1)
      MPLUS=M+1
      IX=0
      IY=N+1
      IF (A(1).GT.A(N)) GOTO 4
    1 MID=(IX+IY)/2
      IF (X.GE.A(MID)) GOTO 2
      IY=MID
      GOTO 3
    2 IX=MID
    3 IF (IY-IX.GT.1) GOTO 1
      GOTO 7
    4 MID=(IX+IY)/2
      IF (X.LE.A(MID)) GOTO 5
      IY=MID
      GOTO 6
    5 IX=MID
    6 IF (IY-IX.GT.1) GOTO 4
    7 NPTS=M+2-MOD(M,2)
      IP=0
      L=0
      GOTO 9
    8 L=-L
      IF (L.GE.0) L=L+1
    9 ISUB=IX+L
      IF ((1.LE.ISUB).AND.(ISUB.LE.N)) GOTO 10
      NPTS=MPLUS
      GOTO 11
   10 IP=IP+1
      T(IP)=A(ISUB)
      D(IP)=F(ISUB)
   11 IF (IP.LT.NPTS) GOTO 8
      EXTRA=NPTS.NE.MPLUS
      DO 14 L=1,M
      IF (.NOT.EXTRA) GOTO 12
      ISUB=MPLUS-L
      D(M+2)=(D(M+2)-D(M))/(T(M+2)-T(ISUB))
   12 I=MPLUS
      DO 13 J=L,M
      ISUB=I-L
      D(I)=(D(I)-D(I-1))/(T(I)-T(ISUB))
      I=I-1
   13 CONTINUE
   14 CONTINUE
      SUM=D(MPLUS)
      IF (EXTRA) SUM=0.5D0*(SUM+D(M+2))
      J=M
      DO 15 L=1,M
      SUM=D(J)+(X-T(J))*SUM
      J=J-1
   15 CONTINUE
      CHUTAB=SUM
      END
CDECK  ID>, CHVDIF.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Bryan Webber
C-----------------------------------------------------------------------
      SUBROUTINE CHVDIF(N,P,Q,R)
C-----------------------------------------------------------------------
C     VECTOR DIFFERENCE (extracted from HERWIG)
C-----------------------------------------------------------------------
      DOUBLE PRECISION P(N),Q(N),R(N)
      INTEGER N,I
      DO 10 I=1,N
   10 R(I)=P(I)-Q(I)
      END
CDECK  ID>, CHVEQU.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Bryan Webber
C-----------------------------------------------------------------------
      SUBROUTINE CHVEQU(N,P,Q)
C-----------------------------------------------------------------------
C     VECTOR EQUALITY (extracted from HERWIG)
C-----------------------------------------------------------------------
      DOUBLE PRECISION P(N),Q(N)
      INTEGER N,I
      DO 10 I=1,N
   10 Q(I)=P(I)
      END
CDECK  ID>, CHVSUM.
*CMZ :-        -17/07/03  18.11.30  by  Peter Richardson
*-- Author :    Bryan Webber
C-----------------------------------------------------------------------
      SUBROUTINE CHVSUM(N,P,Q,R)
C-----------------------------------------------------------------------
C    VECTOR SUM (extracted from HERWIG)
C-----------------------------------------------------------------------
      DOUBLE PRECISION P(N),Q(N),R(N)
      DO 10 I=1,N
   10 R(I)=P(I)+Q(I)
      END
********************************************************************  
****************my** Functions by Sabine
******************** Hossenfelder
******************************************************************

!C#################################################################################
!C######  Used Functions  
!C#################################################################################



!C#################################################################################
!C###### calculate entropy 
!C###### S:=(x,d)->2*Pi*(d+1)/(d+2)*(Mf*RH(x,d))^(d+2) 
!C#################################################################################

 
	function Sf(ddum,Mdum,bhrdum)
      implicit none
      real*8 Sf,Mdum,bhrdum
      integer ddum
      real*8 pi,a0,a1,a2,a3       
      real*8 myM_f,myMrel
      common /constants/ myM_f,myMrel
     
	pi = acos(-1d0)
	a0 = 1d0
	a1 = 0d0
	a2 = 0d0
	a3 = 0d0


      Sf = 2d0*pi*(ddum+1d0)/(ddum+2d0)*
     &(myM_f*bhrdum)**(ddum+2d0)
     &    *(a0 + a1*(myMrel/Mdum) + a2*(myMrel/Mdum)**2d0 
     &    + a3*(myMrel/Mdum)**3d0  ) 
	
      return
      end

 
!C#################################################################################
!C --------  calculate spectral energydensity 
!C---------  n:=(M,omega,d)->exp(S(M-j*omega,d)-S(M,d))*(omega)^(3)*Heaviside(x)
!C#################################################################################


	function Specf(ddum,Mdum,odum,statdum,bhrdum)
      implicit none
      real*8 Mdum,odum,Specf,Sf
      integer ddum,jdum,jend
      real*8 M_f,Mrel,statdum,bhrdum
      common /constants/ M_f,Mrel

!	jend = Int((Mdum-Mrel)/odum)
	
!	if(jend.gt.1d6)then
!	print*,"Do not set omega to zero, otherwise j_max = infinity and I am stuck"
!	endif
!	Specf  = 0d0

!  	do 33 jdum = 1,jend,1
  
	if(odum.lt.(Mdum-Mrel))then
        Specf =   
     &       1/(exp(-Sf(ddum,Mdum-odum,bhrdum)
     &       +Sf(ddum,Mdum,bhrdum))+statdum)
     &       *(odum)**(2d0)
	else
	Specf = 0d0
	endif 
	
!33    continue
      
      return
      end
  




