#define Pi 3.1415926535897932384626433832795029D0

***********************************************************

      Subroutine set_geometry()
      Implicit None

#include "geometry.h"
      
      integer ii_eta,jj_phi
      
      nPhiCEM   = 24
      dPhiCEM   = 2D0*Pi/nPhiCEM
      
      maxEtaCEM = 1.1D0
      nEtaCEM	= 20D0
      dEtaCEM	= 2D0*maxEtaCEM/nEtaCEM

      R_CEM	= 175.26D0

      xTowWidth = 48.48D0

C--- Muon is 2-tower if closer than 1.58 cm from edge of tower
      MuonDZ    = 1.58D0
C    Deta
C    1
C    0	    Dphi -1 0 1
C   -1
C   val of 0 means knock it out
      do ii_eta=0,2
        do jj_phi=0,2
          If(jj_phi.EQ.0 .AND. ii_eta.NE.1) Then
	    knock_region_el(ii_eta,jj_phi)=1
	  Else
	    knock_region_el(ii_eta,jj_phi)=0
          Endif
          If(jj_phi.NE.1) Then
            knock_region_mu(ii_eta,jj_phi)=1
          Else
	    knock_region_mu(ii_eta,jj_phi)=0
	  Endif
	enddo
      enddo
      
      End

***********************************************************

      Subroutine set_Mip()
      Implicit None
      
#include "mip.h"
      
      mip_avg  = -0.533;
      mip_sig  =  0.117;
      mip_zero =  0.01;

      UE_avg   =  0.149;
      
      
      End

***********************************************************
      
      double precision Function getPhi(p)
      Implicit none
      
      double precision p(0:3)
      
      getPhi=atan2(p(2),p(1))
      If(getPhi.LT.0D0) getPhi=getPhi+2D0*Pi
      
      End

***********************************************************
      
      double precision Function getEta(p)
      Implicit none
      
      double precision p(0:3),absp
      
      absp=sqrt(p(1)**2+p(2)**2+p(3)**2)
      If(absp.NE.abs(p(3))) Then
        getEta=1D0/2D0*log((absp+p(3))/(absp-p(3)))
      Else
        getEta=1D20
      Endif
            
      End

***********************************************************
      
      double precision Function getZ(eta,R)
      Implicit none
      
      double precision eta,R
      
      If(eta.LT.20D0) Then
        getZ=R*sinh(eta)
      Else
        getZ=1D20
      Endif
            
      End

***********************************************************

      integer function getiPhi(phi)
      implicit none

#include "geometry.h"
      
      double precision phi,dpphi

      dpphi=phi/dPhiCEM
      getiPhi  = dpphi
      if(getiPhi .GT. nPhiCEM-1)  getiPhi = nPhiCEM-1
      if(getiPhi .LT. 0)          getiPhi = 0

      end

***********************************************************

      integer function getiEta(eta)
      implicit none

#include "geometry.h"
      
      double precision eta

      getieta = -100;

      if(Abs(eta).LT.maxEtaCEM) Then
        getieta=10+eta/dEtaCEM
      endif

      end

***********************************************************

      double precision function minZTow(ieta)
      implicit none
      
#include "geometry.h"      
      
      integer ieta
      double precision mineta
  
      mineta = -maxEtaCEM+ieta*dEtaCEM
      minZTow= R_CEM*sinh(mineta)
  
      end

***********************************************************

      double precision function maxZTow(ieta)
      implicit none
      
#include "geometry.h"      
      
      integer ieta
      double precision maxeta
  
      maxeta = -maxEtaCEM+(ieta+1)*dEtaCEM
      maxZTow= R_CEM*sinh(maxeta)
  
      End

***********************************************************

      double precision function getLocalX(phi)
      implicit none
      
#include "geometry.h"      
      
      double precision phi,phiMod

      phiMod = Abs(mod(phi,dPhiCEM))
      getLocalX   = xTowWidth*(0.5D0 - phiMod/dPhiCEM)

      End

***********************************************************

      double precision function getLocalZ(eta)

      implicit none
      
#include "geometry.h"      
      
      double precision eta,minztower,maxztower,minZTow,maxZTow
      double precision z,getZ,dz1,dz2
      integer ieta,getieta
      
      getLocalZ = -99999.0D0
      ieta = getiEta(eta)

      if(ieta.GE.0 .AND. ieta.LT.20) then
        minztower = minZTow(ieta)
        maxztower = maxZTow(ieta)
        z = getZ(eta , R_CEM)    
        dz1     = z - minztower
        dz2     = maxztower - z
    
C--- dz is positive if the nearest tower edge is to the left
        If(dz1.LT.dz2) then
	  getLocalZ = dz1
        Else
	  getLocalZ = -dz2
	Endif
      endif
  
      end

***********************************************************
C--- returns 0 if the photon in a crack between towers

      double precision function cracks(localx,localz)
      implicit none
      
      double precision localx,localz
      
      cracks = 1.0D0
      if(Abs(localx).GT. 23.1D0) cracks = 0.0D0
      if(Abs(localz).LT.  4.2D0) cracks = 0.0D0

      End
	
***********************************************************

      integer function MuonPhoton(pm, pg)
      implicit none
C--- 0 into the recoil
C--- 1 merge with the muon
C--- 2 knock out
C--- 3 in a crack

#include "geometry.h"
      
      double precision pm,pg      
      double precision meta,mphi,geta,gphi,mu_dz,ph_dz,ph_x
      double precision getEta,getPhi,getLocalZ,getLocalX,cracks
      integer retval,iphi1,iphi2,ieta1,ieta2,dieta,diphi
      integer getiPhi,getiEta

      retval=0
      
      meta=getEta(pm)
      mphi=getPhi(pm)
      geta=getEta(pg)
      gphi=getPhi(pg)

      iphi1=getiPhi(mphi)
      iphi2=getiPhi(gphi)

      ieta1=getiEta(meta)
      ieta2=getiEta(geta)

      mu_dz=getLocalZ(meta)
      ph_dz=getLocalZ(geta)

      ph_x =getLocalX(gphi)
  

      if(ieta1.EQ.-100 .OR. ieta2.EQ.-100) then 
        retval=0
      else 
        if(cracks(ph_x,ph_dz).EQ.0D0) then
	  retval=3
        else
          dieta = ieta1-ieta2
          diphi = iphi1-iphi2
          if (diphi .LT. -12) diphi=diphi+24
          if (diphi .GT.  12) diphi=diphi-24      
C--- merging
C------------------------------------------------------------------
          if(diphi .EQ. 0) then
C--- if em energy is in the same tower as muon, add its energy
	    if(dieta .EQ. 0) then
	      retval = 1  !one-tower muon
	    else if (abs(dieta).EQ.1 .AND. 
     &		     (dieta*mu_dz) .LT. MuonDZ .AND.
     &		     (dieta*mu_dz) .GT. 0D0) then
	      retval = 1  !two-tower muon
            endif
	  endif
C--- Knocking out
C------------------------------------------------------------------      
          if(Abs(diphi).LT.2 .AND.
     &	     Abs(dieta).LT.2 .AND. retval .NE. 1) then 	
	    if(knock_region_mu(1-dieta,1+diphi).EQ.0) retval = 2             
	  endif
        endif
      endif
      
      MuonPhoton=retval
      
      End

*****************************************************************

      integer function ElectronPhoton(pe,pg)
      implicit none
C--- 0 into the recoil
C--- 1 merge with the muon
C--- 2 knock out
C--- 3 in a crack

#include "geometry.h"
      
      double precision pe,pg      
      double precision eeta,ephi,geta,gphi,eledz,phodz,pho_x,lxces
      double precision getEta,getPhi,getLocalZ,getLocalX,cracks
      integer retval,iphi1,iphi2,ieta1,ieta2,dieta,diphi
      integer getiPhi,getiEta

      retval=0
      
      eeta=getEta(pe)
      ephi=getPhi(pe)
      geta=getEta(pg)
      gphi=getPhi(pg)

      iphi1=getiPhi(ephi)
      iphi2=getiPhi(gphi)

      ieta1=getiEta(eeta)
      ieta2=getiEta(geta)

      eledz=getLocalZ(eeta)
      phodz=getLocalZ(geta)

      pho_x =getLocalX(gphi)

      if(ieta1.EQ.-100 .OR. ieta2.EQ.-100) then
        retval=0	
      else    
        if(cracks(pho_x,phodz).EQ.0D0) then
	  retval=3
        else
C--- check the sign vs ces x etc
          dieta=ieta2-ieta1
          diphi=iphi2-iphi1
          if(diphi .LT. -12) diphi=diphi+24
          if(diphi .GT.  12) diphi=diphi-24      
C---merging
C-----------------------------------------------------------------
C--- Algorithm as implemented here :
C--- Merge photons with d-eta = -1,0,+1 in the same phi wedge. 
C--- Exceptions are tower 0 (clusters cannot extend across crack) 
C--- and tower 9 (PEM towers cannot enter CEM clusters).
          if((ieta1.GT.9 .AND. ieta2.GT.9) .OR. 
     &       (ieta1.LT.10 .AND. ieta2.lT.10) )then	
	    if(diphi.EQ.0) then
	      if(dieta.EQ.0) then 
	        retval = 1
	      else if (abs(dieta).EQ.1 .AND.
     &                 (-dieta*eledz).GT.0D0) then
                retval = 1
	      endif
	    endif
	  endif	
C--- Knocking out
C------------------------------------------------------------------    
          if (retval .NE.1) then	
C--- the sign of local x, 
C--- which is falling towards increasing iphi, so we have a -sign
C--- lxces goes from -24 to 24cm when phi increases
	    lxces = -getLocalX(ephi)
	    lxces=lxces/Abs(lxces)	
C--- if dphi>0 means photon on the side of higher phi, 
C--- so electron is closer for positive lxces
C--- if dphi<0 lxces is closer at negative lxces
	    diphi=diphi*lxces	
	    if(Abs(diphi).LT.2 .AND. Abs(dieta).LT.2) then
	      if(knock_region_el(1-dieta,1+diphi).EQ.0) retval=2
            endif
          endif
	endif  
      endif
      ElectronPhoton=retval
      
      End

****************************************************************

C--- Muon energy measured in the EM towers
C--- MIP contribution + underlying event + photons
C--- This routine returns MIP + UE
C--- for photons: check if they need to be added 
C--- using the return value (=1) of the MuonPhoton routine

      double precision function MipE()
      implicit none

#include "mip.h"
      
      double precision ifzero,my_random_number,Gaus_random_number
      
      mipe = UE_avg
  
      ifzero = my_random_number()

      if(ifzero.GT.mip_zero) then !0.01% have zero MIP contribution
        mipe=mipe+10D0**(Gaus_random_number(mip_avg,mip_sig))
      endif	
  
      End

****************************************************************

      integer function MipCutFail(e_em,pt)
      implicit none 
      
      double precision e_em,pt
      integer fail

      fail=0

      if(pt.LE.100) then
        if(e_em.GT.2.0D0) fail=1
      endif
      if(pt.GT.100D0) then
        if(e_em.GT.(2.0D0 + (pt -100.0D0 )*0.0115D0)) fail=1
      endif

      MipCutFail=fail
      
      end

******************************************************************

      double precision function Gaus_random_number(avg,sig)
      implicit none
      
      double precision rd1,rd2,q,avg,sig,my_random_number
      
      q=2D0
      
      do while (q.GT.1)
        rd1=my_random_number()
        rd2=my_random_number()
	q=(2D0*rd1-1D0)**2+(2D0*rd2-1D0)**2
      enddo
      
      q=sqrt(-2D0*log(q)/q)
      
      Gaus_random_number=avg+((2D0*rd1-1D0)*q)*sig
      
      end

******************************************************************

      double precision function my_random_number()
      implicit none
      
C--- specify your favorite random numer generator here
c      double precision my_random_number_generator_function
C      my_random_number=my_random_number_generator_function()
c      Print *,"You have to provide your favourite random"//
c     &        "number generator here"
c      Stop
c   subroutine generates uniformly distributed random no's x(i),i=1,n
      real rand(1)

      call ranmar(rand,1)
      my_random_number=1d0*rand(1)

      end

