 
*************************************************************************
*                                                                       *
*                           S M E A R, 3.02                             *
*                                                                       *
*                  A Program for Detector Simulation                    *
*                  for e+ e- Linear Collider Studies                    *
*                                                                       *
*           Authors: R.Settles, H.Spiesberger and W.Wiedenmann          *
*       Developed from SMEAR,2.0 written by S.Manly for NLC studies     *
*                                                                       *
*              Please report experience, problems, bugs to:             *
*              H.Spiesberger, e-mail: hspiesb@mail.desy.de              *
*                                                                       *
*************************************************************************
 
C...Interface to calling routine
C...call SMEAR for each particle on event record
      SUBROUTINE SMRINT
      IMPLICIT NONE
      INTEGER N,K
      REAL P,V
      COMMON /LUJETS/ N,K(4000,5),P(4000,5),V(4000,5)
      SAVE /LUJETS/
      INTEGER NMXHEP,NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
      DOUBLE PRECISION PHEP,VHEP
      PARAMETER (NMXHEP=2000)
      COMMON /HEPEVT/ NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP)
     &,               JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP)
     &,               PHEP(5,NMXHEP),VHEP(4,NMXHEP)
      SAVE /HEPEVT/
      INTEGER KLU,I,J,CHARGE,EM,EFF,ERR,NCERR
      DATA NCERR /0/
      REAL IV(4),OV(4)
 
      DO 1 I=1,N
      IF (KLU(I,6).EQ.0) THEN
        CHARGE=0
        ELSE
        CHARGE=1
      ENDIF
      IF (IABS(K(I,2)).EQ.11.OR.K(I,2).EQ.22) THEN
          EM=1
        ELSEIF (IABS(K(I,2)).EQ.13) THEN
          EM=-1
        ELSE
          EM=0
      ENDIF
      DO 2 J=1,4
    2 IV(J)=P(I,J)
      CALL SMEAR(CHARGE,EM,IV,OV,EFF,ERR)
      IF (ERR.EQ.1) THEN
        K(I,1)=21
        NCERR=NCERR+1
        IF (NCERR.LT.20) THEN
          WRITE(6,*) 
     & 'SMRINT: Warning: Error in SMEAR, Particle removed'
          WRITE(6,*) ' K(I,2) = ',K(I,2)
          WRITE(6,*) ' IV = ',IV(1),IV(2),IV(3),IV(4)
          WRITE(6,*) ' only 20 first warnings printed '
        ENDIF
        IF (NCERR.EQ.20) THEN
          WRITE(6,*) 
     & 'SMRINT: Warning: Error in SMEAR, Particle removed'
          WRITE(6,*) ' K(I,2) = ',K(I,2)
          WRITE(6,*) ' IV = ',IV(1),IV(2),IV(3),IV(4)
          WRITE(6,*) 'SMRINT: Error has occured 20 times, '
          WRITE(6,*) '       no further warnings printed'
        ENDIF
      ENDIF
      IF (EFF.EQ.0.OR.EFF.EQ.2) K(I,1)=21
      DO 3 J=1,4
    3 P(I,J)=OV(J)
      P(I,5)=0.0
    1 CONTINUE

      RETURN
      END

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

C ****************************************************************
C VERSION 3.0
C ****************************************************************
C Subroutine to smear four momenta of final state charged or 
C neutral particles.
C Smears as a function of theta according to COMMON /SMRCOM/.
C Also sets a flag to allow for different detection efficiencies
C as a function of theta.
C For further details, see comments in block data SMBDDA.
C
C Author: Steve Manly, Yale University
C Written for NLC physics studies, summer 1995
C Modifications by S.Manly, R.Dubois, B.Barakat 
C Version 3.0 by H.Spiesberger (13/09/96)
C
C ****************************************************************
      SUBROUTINE SMEAR(CHARGE,EM,IV,V,EFF,ERR)
      IMPLICIT NONE
C...INPUT to subroutine
      INTEGER CHARGE
C =0 if particle is neutral
C =1 if particle is charged
      INTEGER EM
C =0 if particle is not an electron or photon and energy smear is
C    to be done using hadronic calorimeter resolution
C =1 if particle is an electron or photon and energy smear is to
C    be done using electromagnetic calorimeter resolution
C =-1 if the particle is a muon and momentum smear is to be done
C    using tracking
      REAL IV(4)
C Unsmeared four vector of particle
C IV(1) = INITIAL PX
C IV(2) = INITIAL PY
C IV(3) = INITIAL PZ
C IV(4) = INITIAL E

C...OUTPUT of subroutine
      REAL V(4)
C Smeared four vector of particle
C V(1) = FINAL PX
C V(2) = FINAL PY
C V(3) = FINAL PZ
C V(4) = FINAL E
      INTEGER EFF
C =0 if particle does not pass efficiency cut (i.e., should be 
C       killed by user)
C =1 if particle does pass efficiency cut (i.e., should not be 
C       killed by user)
C =2 if particle is inside beam pipe.
      INTEGER ERR
C =0 if no error occurred during execution
C =1 if error occurred during execution

C...Common block SMRCOM
      REAL AEM,BEM,AHAD,BHAD,PRES,PMSPTS,PMSPHS,PMSPCS,LOWANG,ANGLE
      REAL SIG_HDIR,SIG_EMDIR,CHGEFF,NEUEFF
      INTEGER SEED
      COMMON/SMRCOM/AEM,BEM,AHAD,BHAD,PRES(18),PMSPTS(18),PMSPHS(18),
     $              PMSPCS(18),LOWANG,
     $              ANGLE(19),CHGEFF(18),NEUEFF(18),SIG_HDIR(18),
     $              SIG_EMDIR(18),SEED
      SAVE /SMRCOM/

C...local variables
      REAL PTOT,PTOT2,COSTH,ANG,TEMP,SIGMA,ESMEAR,EPROB,TESTVAL
      REAL SIGMAP,SIGMAH,PSMEAR,PTSMEAR,PHSMEAR,PZSMEAR
      REAL PTPROB,PHPROB,PZPROB,TESTT,TESTH,TESTZ
      REAL IPX,IPY,IPZ,IE,IPT,IPHI,PX,PY,PZ,E,PT,PSMEARD(3),POLD(3)
      REAL PI
      PARAMETER (PI=3.1415927)
      REAL*8 RNDM
      INTEGER BIN,I,LEN
      REAL VCOV,RAN,SMGAUS
      DIMENSION VCOV(2,2)
      INTEGER ISEED(4) 
      DATA ISEED/2789819, 9029301, 54762713, 33498521/

      ERR=0
      IPX=IV(1)
      IPY=IV(2)
      IPZ=IV(3)
      IE=IV(4)
      IF (CHARGE.NE.0.AND.CHARGE.NE.1) THEN
        WRITE(6,*) 'SMEAR: Wrong code for CHARGE: ',CHARGE
        WRITE(6,*) '       Execution stopped'
        STOP
      ENDIF
      IF (EM.NE.0.AND.EM.NE.1.AND.EM.NE.-1) THEN
        WRITE(6,*) 'SMEAR: Wrong code for EM: ',EM
        WRITE(6,*) '       Execution stopped'
        STOP
      ENDIF

C...Find angle for particle to determine proper bin in lookup table
      PTOT2=IPX**2+IPY**2+IPZ**2
      IF (PTOT2.LE.0..OR.IE.LE.0.) THEN
        ERR=1
        GOTO 1000
      ENDIF
      PTOT=SQRT(PTOT2)
      IPT=SQRT(IPX**2+IPY**2)
      COSTH=IPZ/PTOT
      ANG=ACOSD(COSTH)
      IF (ANG.GT.90.) ANG=180.-ANG
      BIN=0
      DO 100 I=1,18
        IF (ABS(ANG).GE.ANGLE(I).AND.ABS(ANG).LT.ANGLE(I+1)) BIN=I
        IF (I.EQ.18) THEN
          IF (ABS(ANG).EQ.90.) BIN=18
        ENDIF
100   CONTINUE		                                       

C...have an error, set flag and return
      IF (BIN.EQ.0) THEN
        ERR=1
        GOTO 1000	
      ENDIF

C...angular acceptance cut about beampipe, set flag and return
      IF (ABS(ANG).LE.LOWANG) THEN
        EFF=2
        PX=0.
        PY=0.
        PZ=0.
        E=0.
        GOTO 1000
      ENDIF

C...Efficiency calculations:
      EFF=1
C...Efficiency calc for neutral particles
      IF (CHARGE.EQ.0) THEN
        TEMP=RNDM(SEED)
        IF (TEMP.GT.NEUEFF(BIN)) EFF=0
C...Efficiency calc for charged particles
      ELSEIF (CHARGE.EQ.1) THEN
        TEMP=RNDM(SEED)
        IF (TEMP.GT.CHGEFF(BIN)) EFF=0
      ENDIF

C...Start energy smear code for photons
      LEN=0
 200  CONTINUE
      IF (CHARGE.EQ.0.AND.EM.EQ.1) THEN
C...find sigma for energy smear distribution
        SIGMA=IE*(AEM/SQRT(IE)+BEM) 
        ESMEAR=SMGAUS(SIGMA,IE,RAN(ISEED(1)),RAN(ISEED(2)))
C...set up output for photons and return
        E=ESMEAR
        PX=IPX*(ESMEAR/IE)
        PY=IPY*(ESMEAR/IE)
        PZ=IPZ*(ESMEAR/IE)
        GOTO 2000
      ENDIF
    
C...start energy smear code for electrons
      IF (CHARGE.EQ.1.AND.EM.EQ.1) THEN
C...find sigma for energy smear distribution
        SIGMA=IE*(AEM/SQRT(IE)+BEM) 
        ESMEAR=SMGAUS(SIGMA,IE,RAN(ISEED(1)),RAN(ISEED(2)))
C...electron direction not smeared, total momentum smeared by 
C   assuming smeared energy and electron mass (in GeV).  
        PSMEAR=ESMEAR**2-(0.000511)**2
        IF (PSMEAR.LE.0.) THEN
          ERR=1
          GOTO 1000
        ENDIF
        PSMEAR=SQRT(PSMEAR)
        PX=(IPX/PTOT)*PSMEAR
        PY=(IPY/PTOT)*PSMEAR
        PZ=(IPZ/PTOT)*PSMEAR
        E=ESMEAR
        GOTO 2000
      ENDIF
    
C...neutral hadronic particle energy smear
      IF (CHARGE.EQ.0.AND.EM.EQ.0) THEN
C...find sigma for energy smear distribution
        SIGMA=IE*(AHAD/SQRT(IE)+BHAD) 
        ESMEAR=SMGAUS(SIGMA,IE,RAN(ISEED(1)),RAN(ISEED(2)))
C...Directional smearing for neutral hadronic particles further down. 
        PSMEAR=ESMEAR**2-(0.135)**2
        IF (PSMEAR.LE.0.) THEN
          ERR=1
          GOTO 1000
        ENDIF
        PSMEAR=SQRT(PSMEAR)
        PX=(IPX/PTOT)*PSMEAR
        PY=(IPY/PTOT)*PSMEAR
        PZ=(IPZ/PTOT)*PSMEAR
        E=ESMEAR
        GOTO 2000
      ENDIF
    
C...Momentum smear for charged particles (other than electrons).
C   New version: smear in pt and phi
      IF ((CHARGE.EQ.1.AND.EM.EQ.0).OR.EM.EQ.-1) THEN
C...find covariance matrix for momentum smear in 1/pt and phi.
        IF (IPT.EQ.0.) GOTO 1001
        VCOV(1,1)=PRES(BIN)**2 
     &    +(PMSPTS(BIN)/IPT/SQRT(SIND(ABS(ANG))))**2
        VCOV(2,2)=PRES(BIN)**2 
     &    +(PMSPHS(BIN)/IPT/SQRT(SIND(ABS(ANG))))**2
        VCOV(1,2)=-(PMSPCS(BIN)/IPT/SQRT(SIND(ABS(ANG))))**2
        VCOV(2,1)=VCOV(1,2)
        IPHI=ACOS(IPX/IPT)
        IF (IPY.LT.0D0) IPHI=2.*PI-IPHI
        CALL SMGAU2(VCOV,1./IPT,IPHI,RAN(ISEED(1)),RAN(ISEED(2)),
     &       RAN(ISEED(3)),RAN(ISEED(4)),PTSMEAR,PHSMEAR)
        PT=1./PTSMEAR
        PX=PT*COS(PHSMEAR)
        PY=PT*SIN(PHSMEAR)
C...z dimension:
210     CONTINUE
        IF (IPZ.EQ.0.) GOTO 1003
        SIGMAP=PMSPTS(BIN)/IPZ/SQRT(SIND(ABS(ANG)))
        PZSMEAR=SMGAUS(SIGMAP,1./IPZ,RAN(ISEED(1)),RAN(ISEED(2)))
        PZ=1./PZSMEAR
C...have smeared momentum components, find smeared energy for
C   these particles by using smeared momentum and assuming
C   mass of particle to be mass of pion in GeV
C...protects from negative and zero momentum errors
        GOTO 1004
1001    PX=0.0
        PY=0.0
        GOTO 210
1003    PZ=0.0
1004    CONTINUE
        PSMEAR=PX**2+PY**2+PZ**2
        IF (PSMEAR.LE.0.) THEN
          ERR=1
          GOTO 1000
        ENDIF
        PSMEAR=SQRT(PSMEAR)
        IF (EM.EQ.-1) THEN
          E=SQRT(PSMEAR**2+0.106**2)
        ELSE
          E=SQRT(PSMEAR**2+0.140**2)
        ENDIF
      ENDIF

C...check negative energy after smearing
2000  CONTINUE
      IF (E.LE.0..AND.EFF.EQ.1) THEN
        LEN=LEN+1
        IF (LEN.LE.100) GOTO 200
        ERR=1
        GOTO 1000
      ENDIF

C...put in directional smearing
      POLD(1)=PX
      POLD(2)=PY
      POLD(3)=PZ
      IF (CHARGE.EQ.0.AND.EM.EQ.0) THEN
C...for neutral hadrons
        CALL SMEARD(POLD,PSMEARD,SIG_HDIR(BIN))
        PX=PSMEARD(1)  
        PY=PSMEARD(2)  
        PZ=PSMEARD(3)  
C...for photons and electrons
      ELSEIF (EM.EQ.1) THEN
        CALL SMEARD(POLD,PSMEARD,SIG_EMDIR(BIN))
        PX=PSMEARD(1)  
        PY=PSMEARD(2)  
        PZ=PSMEARD(3)  
      ENDIF

C...error return and smeared momenta output
1000  CONTINUE
      IF (ERR.NE.0) RETURN
      V(1)=PX
      V(2)=PY
      V(3)=PZ
      V(4)=E

      RETURN
      END

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

C   SUBROUTINE SMEARD(POLD,PSMEARD,SIGMA_DRC)
C-------------------------------------------------------------------------
C   ACTION:
C   This subroutine performs smearing of the direction of momentum 
C   vectors ONLY.
C   Written for NLC physics studies, Winter 1996
C   INPUT: 
C   POLD : (REAL) A three elements 1-D vector representing 
C          the old momentum components x, y and z.
C   SIGMA_DRC: (REAL) width of the angular smearing (in Degrees)
C   OUTPUT:   
C   PSMEARD  : (REAL) A three elements 1-D vector representing
C              the new (smeared) momentum components x, y and z.
C   AUTHOR: 
C   Basem Barakat, Louisian Tech University
C   barakat@phys.latech.edu 318-257-4053
C   METHOD:
C   The overall directional smearing is provided by "swinging"
C   the momentum vector in a cone around the original momentum
C   direction. The distribution of the half angle of the cone
C   GAMMA is assumed to be GAUSSIAN with a width of SIGMA_DRC.
C
C   First we start in a "primed" coordinates system in which the 
C   original momentum vector is completly along the Z_PRM direction. 
C   In this coordinates the unsmeared momentum vector components
C   P_ORG_PRM (PX_ORG_PRM=0., PY_ORG_PRM=0., PZ_ORG_PRM=|P|)
C
C   Then in the new primed coordinates system we smear the direction 
C   of the original momentum vector by an angle GAMMA in a cone
C   around the Z_PRM direction. The resulting smeared momentum
C   components P_SMR_PRM = (PX_SMR_PRM, PY_SMR_PRM, PZ_SMR_PRM)
C
C   Finally, the smeared momentum vector is rotated back to the 
C   original lab coordinate 
C   P_SMR_LAB = (PX_SMR_LAB, PY_SMR_LAB, PZ_SMR_LAB)
C-------------------------------------------------------------------------
      SUBROUTINE SMEARD(POLD,PSMEARD,SIGMA_DRC)
      IMPLICIT NONE
      REAL POLD(3), PSMEARD(3)
      REAL PX_ORG_LAB, PY_ORG_LAB, PZ_ORG_LAB
      REAL PX_ORG_PRM, PY_ORG_PRM, PZ_ORG_PRM
      REAL PX_SMR_PRM, PY_SMR_PRM, PZ_SMR_PRM
      REAL PX_SMR_LAB, PY_SMR_LAB, PZ_SMR_LAB
      REAL THETA_ORG_LAB, PHI_ORG_LAB
      REAL GAMMA, SIGMA_DRC, P_MAG, SMGAUS, P_MAG_NEW
      REAL PHI_SMEAR
      REAL SIGMA_RAD
Chsp
      REAL RAN
Chsp
      INTEGER ISEED(3) 
      DATA ISEED/2789819, 9029301, 54762713/
 
      PX_ORG_LAB=POLD(1)
      PY_ORG_LAB=POLD(2)
      PZ_ORG_LAB=POLD(3)

C...convert into radian
      SIGMA_RAD=SIGMA_DRC/57.29578
 
C...calculate magnitude of momentum vector
      P_MAG=SQRT(PX_ORG_LAB**2+PY_ORG_LAB**2+PZ_ORG_LAB**2)
 
C...calculate polar angle of original momentum vector in lab system
      CALL SMANGL(PX_ORG_LAB,PY_ORG_LAB,PZ_ORG_LAB, 
     +           THETA_ORG_LAB, PHI_ORG_LAB)
 
C...Pick an angle GAMMA (smearing angle) in the primed coordinate system.  
C   PHI_SMEAR is the polar angle PHI in the primed coordinates
      GAMMA=SMGAUS(SIGMA_rad,0.,RAN(ISEED(1)),RAN(ISEED(2)))
      PHI_SMEAR=6.283185307*RAN(ISEED(3))
 
C...Calculate smeared momentum vector componets in the primed coordinates
C   system. Here we assume that the transverse componets x' and y' are
C   equal
      PZ_SMR_PRM=P_MAG*COS(GAMMA)
      PY_SMR_PRM=P_MAG*SIN(GAMMA)*SIN(PHI_SMEAR) 
      PX_SMR_PRM=P_MAG*SIN(GAMMA)*COS(PHI_SMEAR) 
 
C...Rotate the smeared momentum vector back into lab coordinate system. 
      CALL SMROTA(PX_SMR_PRM,PY_SMR_PRM,PZ_SMR_PRM,
     +              PX_SMR_LAB,PY_SMR_LAB,PZ_SMR_LAB,
     +              THETA_ORG_LAB,PHI_ORG_LAB)
      PSMEARD(1)=PX_SMR_LAB
      PSMEARD(2)=PY_SMR_LAB
      PSMEARD(3)=PZ_SMR_LAB
 
      RETURN
      END

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

      SUBROUTINE SMANGL(PX ,PY ,PZ, THETA, PHI)
C...This subroutine calculate the polar angles theta,phi
      IMPLICIT NONE
      REAL PX, PY, PZ, THETA, PHI, PR, COSTH, RP
Chsp
      REAL RAN
Chsp
      INTEGER ISEED
      DATA ISEED/93749813/

      RP=SQRT(PZ**2+PY**2+PX**2)
      COSTH=PZ/RP
      THETA=ACOS(COSTH)
      IF(PY.EQ.0..AND.PX.EQ.0.) THEN
        PHI=6.283185307*RAN(ISEED)
       ELSE
        PHI=ATAN2(PY,PX)
        IF(PHI.LT.0.) PHI=PHI+6.283185307
      ENDIF
      RETURN
      END

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

      FUNCTION SMGAUS(SIGMA,XMEAN,RANDOM1,RANDOM2)
C...By giving two random numbers RANDOM1 and RANDOM2,
C   this function generates a Gaussian distribution
C   with half-width SIGMA around the mean value XMEAN.
      IMPLICIT NONE
      REAL SMGAUS,SIGMA,XMEAN,RANDOM1,RANDOM2
      SMGAUS=1.414213562*SIGMA*SQRT(-ALOG(1.-RANDOM1))
     +  *SIN(6.283185307*RANDOM2)+XMEAN
      END

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

      SUBROUTINE SMGAU2(VCOV,XMEAN1,XMEAN2,R1,R2,R3,R4,X1,X2)
C...By giving four random numbers R1, R2, R3 and R4, this routine 
C   generates two variables X1, X2 which are Gaussian distributed 
C   with covariance matrix VCOV around the mean values XMEAN1, XMEAN2.
      IMPLICIT NONE
      REAL VCOV,XMEAN1,XMEAN2,R1,R2,R3,R4,X1,X2
      REAL G1,G2,CCOV
      DIMENSION VCOV(2,2),CCOV(2,2)
      G1=1.414213562*SQRT(-ALOG(1.-R1))*SIN(6.283185307*R2)
      G2=1.414213562*SQRT(-ALOG(1.-R3))*SIN(6.283185307*R4)
      CCOV(1,1)=SQRT(VCOV(1,1))
      CCOV(2,1)=VCOV(2,1)/CCOV(1,1)
      CCOV(2,2)=SQRT(ABS(VCOV(2,2)-VCOV(2,1)**2/VCOV(1,1)))
      X1=G1*CCOV(1,1)+XMEAN1
      X2=G1*CCOV(2,1)+G2*CCOV(2,2)+XMEAN2
      RETURN
      END

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

      SUBROUTINE SMROTA(PXO,PYO,PZO,PXN,PYN,PZN,THETA,PHI)
C...This subroutine calculates componets of a vector in a rotated 
C   coordinate system
C   Old componets:  PXO,PYO,PZO
C   New componets:  PXN,PYN,PZN
C   Rotation angles: THETA, PHI
      IMPLICIT NONE
      REAL PXO,PYO,PZO,PXN,PYN,PZN,THETA,PHI
      REAL SIPHI, COPHI, SITH, COTH
      SIPHI=SIN(PHI)
      COPHI=COS(PHI)
      SITH=SIN(THETA)
      COTH=COS(THETA)
      PXN=-PXO*SIPHI-PYO*COTH*COPHI+PZO*SITH*COPHI
      PYN=PXO*COPHI-PYO*COTH*SIPHI+PZO*SITH*SIPHI
      PZN=PYO*SITH+PZO*COTH
      RETURN
      END
 
*************************************************************************

      BLOCK DATA SMBDDA
      IMPLICIT NONE
C
C User will want to adjust the quantities in COMMON /SMRCOM/ as desired.
C
C Smear energy of photons and electrons as Gaussian with half-width,
C sigma, given by the following expression:
C sigma/IE = AEM/sqrt(IE) + BEM 
C 
C Direction of photons and electrons smeared in cone about the original 
C direction with the half angle of the cone Gaussian distributed with
C half width SIG_EMDIR (corresponds to 10mrad).
C For electrons, the magnitude of the output momentum is determined
C from the smeared calorimeter energy assuming the mass of the particle
C is that of an electron.  
C
C Smear energy of neutral particles interacting hadronically as 
C Gaussian with half-width, sigma, given by the following expression:
C sigma/IE = AHAD/sqrt(IE) + BHAD 
C
C Neutral hadronic particle total momentum smeared by assuming smeared
C energy and pizero mass (in GeV). The momentum components are smeared
C accordingly. Direction smeared in cone about the original direction with
C the half angle of the cone Gaussian distributed with half width SIG_HDIR
C (corresponding to 15mrad).
C
C For charged particles (other than electrons) the smeared energy
C is determined from the momentum assuming the particle is a pion.
C
C Momentum components of charged particles not interacting electro-
C magentically are smeared as follows:
C 1. Transverse momentum Pt and azimuth Phi are smeared together, 
C taking into account their correlation (see Gluckstern, Nucl.Inst.Meth.
C 24 (1963) 381). The components of the covariance matrix VCOV(i,j) are 
C parametrized as follows:
C for 1/Pt:  VCOV(1,1) = PRES**2 + (PMSPTS/(Pt*sqrt(sin(ANG)))**2
C for Phi:   VCOV(2,2) = PRES**2 + (PMSPHS/(Pt*sqrt(sin(ANG)))**2
C for the correlation: 
C            VCOV(1,2) = PRES**2 + (PMSPCS/(Pt*sqrt(sin(ANG)))**2
C i.e. 1/Pt and Phi are Gaussian of half-width given by sqrt(VCOV(1,1)),
C sqrt(VCOV(2,2)), resp. and correlation VCOV(1,2). PRES and PMSPxS are
C given as function of theta (one value given for each 5 degree increment
C in the detector). ANG is the angle.
C 2. The z-component is smeared by projecting the transverse momentum
C i.e. according to a Gaussian in 1/Pz with half-width
C Pzerror/Pz**2 = PMSPTS/(Pz*sqrt(sin(ANG))).
C
C Angular dependence is put into this version in 5 degree increments
C (with respect to the beampipe). There is symmetry about the midplane.
C ANGLE is the array used for angular dependence bookkeeping in SMEAR.
C User should not change this.
C
C Acceptance hole around beampipe determined by value of LOWANG.
C Tracks with an angle with respect to the beampipe that is less that 
C LOWANG will have EFF=2. Output 4-vector will be zeroed out.
C
C No attempt is made in this version to take into account energy 
C cluster overlaps for closely spaced particles (i.e., those
C in narrow jets)
C
C Charged particle efficiencies CHGEFF are given as a function of angle 
C (in 5 degree increments) with respect to the beampipe. Symmetry with 
C respect to the detector midplane is assumed.
C
C Neutral particle efficiencies NEUEFF are given as a function of angle 
C (in 5 degree increments) with respect to the beampipe. Symmetry with 
C respect to the detector midplane is assumed.
C
C SMRCOM written to be used with SUBROUTINE SMEAR, which is a program
C designed to smear particle 4 momenta for NLC physics studies
C 7/27/95
C Steve Manly - in%"manly@yalph2.physics.yale.edu"
C Yale University
C
C Major updates by H.Spiesberger - e-mail: hspiesb@mppmu.mpg.de
C Max-Planck-Institut Muenchen
C 13/09/96

C...Common block SMRCOM
      REAL AEM,BEM,AHAD,BHAD,PRES,PMSPTS,PMSPHS,PMSPCS,LOWANG,ANGLE
      REAL SIG_HDIR,SIG_EMDIR,CHGEFF,NEUEFF
      INTEGER SEED
      COMMON/SMRCOM/AEM,BEM,AHAD,BHAD,PRES(18),PMSPTS(18),PMSPHS(18),
     $              PMSPCS(18),LOWANG,
     $              ANGLE(19),CHGEFF(18),NEUEFF(18),SIG_HDIR(18),
     $              SIG_EMDIR(18),SEED
      SAVE /SMRCOM/

      DATA SEED/9454837/
C Seed is the random number seed used in the smearing

C original data:
C      DATA AEM/.12/
C      DATA BEM/.01/
C 1TeV detector:
      DATA AEM/.10/
      DATA BEM/.01/

C original data:
C      DATA AHAD/.45/
C      DATA BHAD/.02/
C 1TeV detector:
      DATA AHAD/.50/
      DATA BHAD/.04/

      DATA ANGLE/0.,5.,10.,15.,20.,25.,30.,35.,40.,45.,
     $           50.,55.,60.,65.,70.,75.,80.,85.,90./

C original contents
C      DATA PRES/.0002,.0002,.0002,.0002,.0002,.0002,.0002,.0002,.0002,
C     $          .0002,.0002,.0002,.0002,.0002,.0002,.0002,.0002,.0002/

C for 1 TeV detector (R.S.)
      DATA PRES/312.,.045, .0078,.0027,.0013,.00067,.00039,.00024,
     $         .0002,.0002,.0002,.0002,.0002,.0002, .0002, .0002,
     $         .0002,.0002/

C...Momentum resolution due to multiple scattering
C...for 1/pt 
      DATA PMSPTS/.0015,.0015,.0015,.0015,.0015,.0015,.0015,
     $            .0015,.0015,.0015,.0015,.0015,.0015,.0015,.0015,
     $            .0015,.0015,.0015/
C...for phi
      DATA PMSPHS/.00044,.00044,.00044,.00044,.00044,.00044,.00044,
     $            .00044,.00044,.00044,.00044,.00044,.00044,.00044,
     $            .00044,.00044,.00044,.00044/
C...for 1/pt - phi correlation
      DATA PMSPCS/.00030,.00030,.00030,.00030,.00030,.00030,.00030,
     $            .00030,.00030,.00030,.00030,.00030,.00030,.00030,
     $            .00030,.00030,.00030,.00030/

C original data:
C      DATA LOWANG/15./
C 1TeV detector:
      DATA LOWANG/8.1/

C...Efficiencies
      DATA CHGEFF/.98,.98,.98,.98,.98,.98,.98,.98,.98,.98,.98,.98,
     $            .98,.98,.98,.98,.98,.98/
      DATA NEUEFF/.98,.98,.98,.98,.98,.98,.98,.98,.98,.98,.98,.98,
     $            .98,.98,.98,.98,.98,.98/
      DATA SIG_HDIR/0.86,0.86,0.86,0.86,0.86,0.86,0.86,0.86,0.86,
     $              0.86,0.86,0.86,0.86,0.86,0.86,0.86,0.86,0.86/
      DATA SIG_EMDIR/0.57,0.57,0.57,0.57,0.57,0.57,0.57,0.57,0.57,
     $               0.57,0.57,0.57,0.57,0.57,0.57,0.57,0.57,0.57/
      END

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

C...Reset detector smearing parameters for LEP/SLD detector
      SUBROUTINE SMBDD2
      IMPLICIT NONE
C...Common block SMRCOM
      REAL AEM,BEM,AHAD,BHAD,PRES,PMSPTS,PMSPHS,PMSPCS,LOWANG,ANGLE
      REAL SIG_HDIR,SIG_EMDIR,CHGEFF,NEUEFF
      INTEGER SEED
      COMMON/SMRCOM/AEM,BEM,AHAD,BHAD,PRES(18),PMSPTS(18),PMSPHS(18),
     $              PMSPCS(18),LOWANG,
     $              ANGLE(19),CHGEFF(18),NEUEFF(18),SIG_HDIR(18),
     $              SIG_EMDIR(18),SEED
      SAVE /SMRCOM/
      INTEGER I

C LEP/SLD detector:
      AEM=.20
      BEM=.01
      AHAD=.90
      BHAD=.02
      PRES(1)=318.
      PRES(2)=.102
      PRES(3)=.020
      PRES(4)=.0075
      PRES(5)=.0036
      PRES(6)=.0019
      PRES(7)=.0011
      PRES(8)=.0007
      PRES(9)=.0006
      PRES(10)=.0006
      PRES(11)=.0006
      PRES(12)=.0006
      PRES(13)=.0006
      PRES(14)=.0006
      PRES(15)=.0006
      PRES(16)=.0006
      PRES(17)=.0006
      PRES(18)=.0006
      DO 1 I=1,18
      PMSPTS(I)=0.0050
      PMSPHS(I)=0.00061
    1 PMSPCS(I)=0.00053
      LOWANG=15.
      RETURN
      END

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

      FUNCTION RAN()
      IMPLICIT NONE
      REAL RAN
      REAL*8 RNDM
      RAN=RNDM()
      RETURN
      END
