C   20/01/83 705071225  MEMBER NAME  TRLGS6   (S)           FORTRAN
      SUBROUTINE TRLGS6(PV,R,*)
C
C      THIS ROUTINE TRACKS PARTICLES THROUGH THE LEAD GLASS
C      ARRAYS AND ACCUMULATES PULSEHEIGHTS
C      PV = INCOMING 4 VECTOR IS CHANGED DUE TO ENERGY LOSS
C      R  = 3 VECTOR POINTING TO THE FACE OF THE LEAD GLASS
C      RETURN 1 IS USED IF THE INCIDENT MOMENTUM IS LESS THAN
C               10 MEV
C
C      CHANGE       04/02/83    K.MEIER
C      CHANGE OF EXPECM TO INCORPORATE THE SF6 BLOCKS
C                   25/10/85    N.MAGNUSSEN
C      THIS IS ALMOST IDENTICAL TO TRLG1 BUT THE 3-DIMENSIONAL
C      SHOWER INTEGRATION IS DONE ONLY FOR ELECTRONS AND GAMMAS
C      WITH E > 100 MEV (STANDARD PROCEDURE IS USED FOR LOW ENERGIES)
C      CHANGE       20/09/83    K.MEIER
C      ERROR IN CHARGE CALCULATION CORRECTED
C    PRINT STATEMENTS COMMENTED OUT  11.6.1986    J.OLSSON
C    ERROR IN RATION ESEEN/ENERGY CORRECTED 3.7.86  19.30  J.OLSSON
C    ERROR IN NOT UPDATING R-VECTOR FOR ENDCAP PHOTONS, LEADING TO
C    DOUBLE TRACKING OF ENDCAP PHOTONS AND ELECTRONS, CORRECTED
C    LIKEWISE THE SCALING OF BLOCK ENERGIES DEPENDENT ON READOUT
C    THRESHOLD REMOVED, SINCE THIS IS HANDLED CORRECTLY IN RDALGN
C                                  20.08.1986     J.OLSSON
C    UPDATE ENDCAP POS.VECTOR WITH 20 % EXTRA, FOR PRECISION IN TRLGL
C                                   4.9.86   J.OLSSON
C
C    CALL TO VEGAS REPLACED BY CALL TO WEGAS (COMMON TO TRLGS6,TRLGSH)
C    FNORM FACTORS UPDATED ACC. TO F11MAG.CLUS VALUES
C    RENAME EXPECM TO EXPEC6, TO ENABLE RUNNING AS TRLGS6 TOGETHER WITH
C    TRLGSH, WHICH CALLS ANOTHER EXPECM...    17.11.1986  J.O.
C      CHANGES      19/01/87    J.OLSSON/N.MAGNUSSEN
C      2. FNORM IN EXPECM CHANGED TO 59000 AND 21220 (INTEGRATION VALUE)
C
C   CHANGE OF CODE IN SF65, TO PREVENT NEG. ARGUMENT IN ALOG
C                                       30.1.87    N.M / J.O
C   CHANGE OF CODE IN SF56, FOR SIMILAR REASONS..
C                                        2.2.87    N.M / J.O
C   CHANGE OF DIMENSION IN LISTBL,BLCFRC
C                                        1.4.87    D.P / J.O
C   AVGSUM NORM IN ENERGY SUM
C                                       13.4.87    D.P / J.O
C   PREVENT DIVIDE CHECK FOR AVGSUM=0.
C                                        2.5.87          J.O
C
C
      COMMON/CGEO1/BKGAUS, RPIP,DRPIP,XRLPIP, RBPC,DRBPC,XRLBPC,
     *             RITNK,DRITNK,XRLTKI, R0ROH,DR0ROH,XR0ROH,
     *             R1ROH,DR1ROH,XR1ROH, R2ROH,DR2ROH,XR2ROH,
     *             R3ROH,DR3ROH,XR3ROH, ROTNK,DROTNK,XRLTKO,
     *             RTOF,DRTOF,XRTOF, RCOIL, DRCOIL, XRCOIL,
     *             ZJM,DZJM,XRZJM, ZJP,DZJP,XRZJP,
     *             ZTKM,DZTKM,XRZTKM, ZTKP,DZTKP,XRZTKP,
     *             ZBPPL,ZBPMI,ZTOFPL,ZTOFMI,
     *             XRJETC,
     *             RLG,ZLGPL,ZLGMI,OUTR2,CTLIMP,CTLIMM,DELFI,
     *             BLXY,BLZ,BLDEP,ZENDPL,ZENDMI,DEPEND,
     *             XHOL1,XHOL2,YHOL1,YHOL2
      COMMON/CLGHIT/AMPLG(3000)
      DIMENSION RR(3),LISTBL(100),BLCFRC(100),VAR(2),DC(3)
      LOGICAL * 1 LFLAG
      COMMON/CFLAG/LFLAG(10)
C        LFLAG(1) = SMEAR GAMMA AND ELECTRON ENERGIES
C        LFLAG(2) = GAMMA CONVERSION IN OUTER TANK AND COIL (TRKGAM)
C        LFLAG(3) = ABSORPTION LOSSES
C        LFLAG(4) = 3 DIM SHOWER PROFILE FIT TO EGS CODE SLOW!!!!!
C        LFLAG(5) = .TRUE.   -->  WITH VERTEX CHAMBER TRACKING
C                 = .FALSE.  -->  WITHOUT VERTEX CHAMBER TRACKING
C                                   BUT OLD BEAM PIPE GEOMETRY AND
C                                 BEAM PIPE COUNTERS (BEFORE MAI 84)
C
C                     ===> DEFAULT .TRUE.,.FALSE.,.TRUE.,.FALSE.,.FALSE.
C
      DATA TWOPI/6.28319/
      REAL RADL/23.9/
*** PMF 15/10/99      DIMENSION R(1),PV(1)
      DIMENSION R(*),PV(*)
*** PMF (end)
      DIMENSION R0(3)
      DATA ICALL /0/
C
      IF(ICALL.EQ.1) GOTO 1111
      WRITE(6,1113)
1113  FORMAT(' ===>> TRLGS6: DOUBLE TRACKING CORRECTED VERSION 4.9.86 '/
     $ '                    REVISED VERSION FROM 19.01.1987...')
      ICALL = 1
C
C      FIND INDEX FOR ENERGY LOSS
C       LOSS = 0  PHOTON
C       LOSS = 1  ELECTRON
C       LOSS = 2  MINIMUM IONIZING PARTICLE
C
 1111 IF(PV(6).LT.0.01) RETURN 1
C
      OUTR=OUTR2*OUTR2
      ELOSS=0.
      ELS  =0.
      EDEPOT = 0.
      LOSS=2
      IF(PV(5).LT.1.E-03) LOSS=1
      IF(PV(5).EQ.0.)     LOSS=0
C
C                                           ENERGIES IN MEV
      ENERGY=PV(4)*1000.
      EINC = ENERGY / 1000.
C
C                                           THE DIRECTION OF THE TRACKIN
C                                           VECTOR IS NOT CHANGED
C                                           WHEN TRACKED THROUGH THE
C                                           BLOCKS
      PV(6)=SQRT(PV(1)**2+PV(2)**2+PV(3)**2)
      PVX=PV(1)/PV(6)
      PVY=PV(2)/PV(6)
      PVZ=PV(3)/PV(6)
      PVXY=PVX**2+PVY**2
      IF(PVXY.LT.1.E-5) RETURN
      DC(1) = PVX
      DC(2) = PVY
      DC(3) = PVZ
      X1=0.
      R1=R(1)**2+R(2)**2
      RADIUS = SQRT( R1 )
      IF( .NOT.
     *  ( ZLGMI+5..LE.R(3).AND.R(3).LE.ZLGPL-5..AND.RADIUS.GT.RLG-5. ))
     *                             GO TO 1000
C
C                                           CENTRAL PART
C
C
C                                           ANTI ENERGY CORRECTION
C                                           ABSORPTION IN COIL
      IBE = 0
      IF( LOSS .LT.2 .AND. LFLAG(3) ) ENERGY=ELGOBS(IBE,DC,EINC)*1000.
C                                           ENERGY SMEARING
      ENERGY = AMAX1( ENERGY, 1. )
      EINC = ENERGY/1000.
      IF( LOSS .LT.2 .AND. LFLAG(1) ) ENERGY=ELSMR(IBE,DC,EINC)*1000.
C
C                                           CALCULATE THICKNESS OF SLAB
C                                           OF LEAD GLASS CORRESPONDING
C                                           TO 1CM TRACK LENGTH
C
      PVR=R(1)*PVX+R(2)*PVY
      DELR=PVR/RADIUS
      TRKL=BLDEP/DELR
C
C                                           DETERMINE STEP SIZE
      NSTEP=BLDEP/(10.*DELR)+1
      XL=10.
      DX=XL/RADL
C
C                                           ENERGY LOSS FOR HADRONS
      IF( LOSS.LT.2.AND.ENERGY.GT.100.) GOTO 500
      IF( LOSS.LT.2 ) GO TO 26
C
C     CALCULATE EXIT POINT AT THE BACK SURFACE OF THE LEAD GLASS
C
      R0(1)=R(1)+PVX*TRKL
      R0(2)=R(2)+PVY*TRKL
      R0(3)=R(3)+PVZ*TRKL
C
C     CALCULATE PATH LENGTH THROUGH LIGHT GUIDE (IFLG = 2)
C
      IFLG = 2
      CALL PLLIGD(PLLI,R(1),R(2),R(3),R0(1),R0(2),R0(3),IFLG)
      ICHG=IFIX(PV(7))
      ITYP=IFIX(PV(8))
      CALL HDECLS(EDEPOT,ICHG,ITYP,PV,IBE,TRKL,PLLI)
      ED=EDEPOT
      EDEPOT=10.*EDEPOT/TRKL*1000.
C     WRITE(6,600) ICHG,ITYP,PV(4),IBE,TRKL,PLLI,ED,EDEPOT
  600 FORMAT(1H ,2I3,F8.3,I5,2F8.2,2X,2F8.3)
   26 CONTINUE
      N0 = 1
C
C                                           DEPTH OF CONVERSION
C                                           FOR PHOTONS IN THE LEAD GLAS
C                                           ASSUME UNIFORM CONVERSION
C                                           PROBABILITY WITHIN
C                                           THE FIRST 2 X0
      IF(LOSS.NE.0) GO TO 23
      X0=51.*RN(DUM)
      N0=X0/(10.*DELR)+1
   23 CONTINUE
C
      ISTEP=0
   21 ISTEP=ISTEP+1
C
C                                           CALCULATE INDEX FOR BLOCKS
C                                           WHICH ARE HIT
      PHI=ATAN2(R(2),R(1))
      IF(PHI.LT.0.) PHI=TWOPI+PHI
      NFI=PHI/DELFI
      IF(ABS(R(3)).GT.ZLGPL) GO TO 22
      DZ=R(3)-ZLGMI
      NZ=DZ/BLZ+1
      NBL=NZ+32*NFI
C
C                                           ACCUMULATE PULSEHEIGHTS
C                                           IN BLOCKS
      IF(LOSS.GE.2) GO TO 25
C                                           PHOTONS AND E
      EDEPOT=0.
      IF(ISTEP.LT.N0) GO TO 24
      X1=X1+DX
      EDEPOT=SHOWR(ENERGY,X1,DX)
      IF(X1.GT.9.) EDEPOT=ENERGY-ELOSS
      ELOSS=ELOSS+EDEPOT
      IF(ELOSS.GT.ENERGY) GO TO 24
      AMPLG(NBL)=AMPLG(NBL)+EDEPOT
      IF(X1.GT.9.) RETURN
      GO TO 24
C
C                                           HADRONS
   25 AMPLG(NBL)=AMPLG(NBL)+EDEPOT
      ELS=ELS+EDEPOT
      IF(ELS.GE.ENERGY) GO TO 22
C
C                                           PROPAGATE TRACK
   24 R(1)=R(1)+XL*PVX
      R(2)=R(2)+XL*PVY
      R(3)=R(3)+XL*PVZ
C
C                                           CHECK WHETHER PARTICLE IS
C                                           STILL IN THE LEAD GLASS
      R1=R(1)**2+R(2)**2
      IF(R1.GT.OUTR) GO TO 22
      IF(ZLGPL.LT.R(3).OR.R(3).LT.ZLGMI) GO TO 22
      GO TO 21
   22 CONTINUE
      IF(ELOSS.EQ.0.) RETURN
      RX=(PV(4)-ELOSS/2000.)/PV(4)
      IF(ELOSS.GE.ENERGY) RX=0.001
      IF(ELS  .GE.ENERGY) RX=0.001
      RX=ABS(RX)
      DO 1004 I4=1,4
 1004 PV(I4)=PV(I4)*RX
      PV(6)=PV(6)*RX
      RETURN
C
C                                           END CAP LEAD GLASS
C                                           FOLLOW PHOTON TO THE FACE
C                                           OF THE END CAP COUNTERS
 1000 CONTINUE
      RABS=0.
      DO 1001 I3=1,3
 1001 RABS=RABS+R(I3)**2
      RABS=SQRT(RABS)
      DO 1002 I3=1,3
 1002 R(I3)=R(I3)/RABS
      DISTZ=ZENDPL
      IF(R(3).LT.0) DISTZ=ABS(ZENDMI)
      ZLAM=ABS(DISTZ/R(3))
      DO 1003 I3=1,3
 1003 R(I3)=ZLAM*R(I3)
      CALL ENDCLG(R,NBL,*9000)
C
C                                           ANTI ENERGY CORRECTION
C                                           ABSORPTION IN TANK ENDPLATE
      IBE = 1
      IF( R(3) .LT. 0. ) IBE = -1
      IF( LOSS .LT.2 .AND. LFLAG(3) ) ENERGY=ELGOBS(IBE,DC,EINC)*1000.
C                                           ENERGY SMEARING
      ENERGY = AMAX1( ENERGY, 1. )
      EINC = ENERGY/1000.
      IF( LOSS .LT.2 .AND. LFLAG(1) ) ENERGY=ELSMR(IBE,DC,EINC)*1000.
C
C                                           CALCULATE DELZ OF SLAB
C                                           OF LEAD GLASS CORRESPONDING
C                                           1CM TRACK LENGTH IN LEAD GLA
      DELZ=ABS(PVZ*10.)
      TRKL=DEPEND/ABS(PVZ)
C
C                                           DETERMINE STEP SIZE
      NSTEP=DEPEND/DELZ +1
      XL=10.
      DX=XL/RADL
C
C                                           ENERGY LOSS FOR HADRONS
      IF( LOSS.LT.2.AND.ENERGY.GT.100.) GOTO 700
      IF( LOSS.LT.2 ) GO TO 260
      PLLI=0.
      ICHG=IFIX(PV(7))
      ITYP=IFIX(PV(8))
      CALL HDECLS(EDEPOT,ICHG,ITYP,PV,IBE,TRKL,PLLI)
      ED=EDEPOT
      EDEPOT=10.*EDEPOT/TRKL*1000.
C     WRITE(6,600) ICHG,ITYP,PV(4),IBE,TRKL,PLLI,ED,EDEPOT
  260 CONTINUE
      N0 = 1
C
      IF(LOSS.NE.0) GO TO 230
C
C                                           DEPTH OF CONVERSION
C                                           FOR PHOTONS IN THE LEAD GLAS
C                                           ASSUME UNIFORM CONVERSION
C                                           PROBABILITY WITHIN
C                                           THE FIRST 2 X0
      IF(LOSS.NE.0) GO TO 230
      X0=51.*RN(DUM)
      N0=X0/DELZ +1
  230 CONTINUE
C
      ISTEP=0
  210 ISTEP=ISTEP+1
C
C                                           CALCULATE INDEX FOR
C                                           BLOCKS WHICH ARE HIT AND
C                                           CHECK WHETHER PARTICLE IS
C                                           STILL IN THE LEAD GLASS
      CALL ENDCLG(R,NBL,*9000)
C
C                                           ACCUMULATE PULSEHEIGHTS
      IF(LOSS.GE.2) GO TO 250
C                                           PHOTONS AND E
      EDEPOT=0.
      IF(ISTEP.LT.N0) GO TO 240
      X1=X1+DX
      EDEPOT=SHOWR(ENERGY,X1,DX)
      IF(X1.GT.8.) EDEPOT=ENERGY-ELOSS
      ELOSS=ELOSS+EDEPOT
      IF(ELOSS.GT.ENERGY) GO TO 240
      AMPLG(NBL)=AMPLG(NBL)+EDEPOT
      IF(X1.GT.8.) RETURN
      GO TO 240
C                                           HADRONS
  250 AMPLG(NBL)=AMPLG(NBL)+EDEPOT
      ELS=ELS+EDEPOT
      IF(ELS.GE.ENERGY) GO TO 9000
C                                           PROPAGATE TRACK
  240 R(1)=R(1)+XL*PVX
      R(2)=R(2)+XL*PVY
      R(3)=R(3)+XL*PVZ
      GO TO 210
C
C                                           DE/DX ENERGY LOSS
 9000 IF(ELOSS.EQ.0.) RETURN
      RX=(PV(4)-ELOSS/2000.)/PV(4)
      IF(ELOSS.GE.ENERGY) RX=0.001
      IF(ELS  .GE.ENERGY) RX=0.001
      RX=ABS(RX)
      DO 9001 I4=1,4
 9001 PV(I4)=PV(I4)*RX
      PV(6)=PV(6)*RX
C
      RETURN
C
C  500   CENTRAL PART ELECTRONS AND PHOTONS
C  700   ENDCAP       ELECTRONS AND PHOTONS
C
C      CENTRAL PART
C      CALCULATE BLOCKS WHICH MAY SEE PART OF THE EL. MG. SHOWER
C
  500 XL = 40.
      NO = 1
      RR(1) = R(1)
      RR(2) = R(2)
      RR(3) = R(3) + 100.
C
C      CHECK WHETHER PARTICLE IS STILL IN THE LEAD GLASS
C
      R1=RR(1)**2+RR(2)**2
      IF(R1.GT.OUTR) GO TO 510
      IF(ZLGPL.LT.RR(3).OR.RR(3).LT.ZLGMI) GO TO 510
  501 PHI=ATAN2(RR(2),RR(1))
      IF(PHI.LT.0.) PHI=TWOPI+PHI
      NFI=PHI/DELFI
      IF(ABS(RR(3)).GT.ZLGPL) GO TO 510
      DZ=RR(3)-ZLGMI
      NZ=DZ/BLZ+1
      NBL=NZ+32*NFI
C
C     CHECK WHETHER NBL IS ALREADY IN THE BLOCK LIST IF NOT
C     ADD THE BLOCK NUMBER TO THE LIST
C
      IF(NO.EQ.1) GO TO 505
      DO 504 N=1,NO
      IF(LISTBL(N).EQ.NBL) GO TO 503
  504 CONTINUE
      NO = NO + 1
  505 LISTBL(NO)=NBL
      NPL = NFI+1
      NMI = NFI-1
      IF(NFI.EQ. 0) NMI = 83
      IF(NFI.EQ.83) NPL =  0
      LISTBL(NO+1) =NZ+32*NPL
      LISTBL(NO+2) =NZ+32*NMI
      NO = NO + 2
  503 RR(1)=RR(1)+XL*PVX
      RR(2)=RR(2)+XL*PVY
      RR(3)=RR(3)+XL*PVZ
C
C      CHECK WHETHER PARTICLE IS STILL IN THE LEAD GLASS
C
      R1=RR(1)**2+RR(2)**2
      IF(R1.GT.OUTR) GO TO 510
      IF(ZLGPL.LT.RR(3).OR.RR(3).LT.ZLGMI) GO TO 510
      GO TO 501
C
  510 RR(1) = R(1)
      RR(2) = R(2)
      RR(3) = R(3) - 100.
C
C      CHECK WHETHER PARTICLE IS STILL IN THE LEAD GLASS
C
      R1=RR(1)**2+RR(2)**2
      IF(R1.GT.OUTR) GO TO 520
      IF(ZLGPL.LT.RR(3).OR.RR(3).LT.ZLGMI) GO TO 520
C
  511 PHI=ATAN2(RR(2),RR(1))
      IF(PHI.LT.0.) PHI=TWOPI+PHI
      NFI=PHI/DELFI
      IF(ABS(RR(3)).GT.ZLGPL) GO TO 520
      DZ=RR(3)-ZLGMI
      NZ=DZ/BLZ+1
      NBL=NZ+32*NFI
C
C     CHECK WHETHER NBL IS ALREADY IN THE BLOCK LIST IF NOT
C     ADD THE BLOCK NUMBER TO THE LIST
C
      IF(NO.EQ.1) GO TO 515
      DO 514 N=1,NO
      IF(LISTBL(N).EQ.NBL) GO TO 513
  514 CONTINUE
      NO = NO + 1
  515 LISTBL(NO)=NBL
      NPL = NFI+1
      NMI = NFI-1
      IF(NFI.EQ. 0) NMI = 83
      IF(NFI.EQ.83) NPL =  0
      LISTBL(NO+1) =NZ+32*NPL
      LISTBL(NO+2) =NZ+32*NMI
      NO = NO + 2
  513 RR(1)=RR(1)+XL*PVX
      RR(2)=RR(2)+XL*PVY
      RR(3)=RR(3)+XL*PVZ
C
C      CHECK WHETHER PARTICLE IS STILL IN THE LEAD GLASS
C
      R1=RR(1)**2+RR(2)**2
      IF(R1.GT.OUTR) GO TO 520
      IF(ZLGPL.LT.RR(3).OR.RR(3).LT.ZLGMI) GO TO 520
C
      GO TO 511
C
  520 RR(1) = R(1)
      RR(2) = R(2)
      RR(3) = R(3)
C
C      CHECK WHETHER PARTICLE IS STILL IN THE LEAD GLASS
C
      R1=RR(1)**2+RR(2)**2
      IF(R1.GT.OUTR) GO TO 530
      IF(ZLGPL.LT.RR(3).OR.RR(3).LT.ZLGMI) GO TO 530
C
  521 PHI=ATAN2(RR(2),RR(1))
      IF(PHI.LT.0.) PHI=TWOPI+PHI
      NFI=PHI/DELFI
      IF(ABS(RR(3)).GT.ZLGPL) GO TO 530
      DZ=RR(3)-ZLGMI
      NZ=DZ/BLZ+1
      NBL=NZ+32*NFI
C
C     CHECK WHETHER NBL IS ALREADY IN THE BLOCK LIST IF NOT
C     ADD THE BLOCK NUMBER TO THE LIST
C
      IF(NO.EQ.1) GO TO 525
      DO 524 N=1,NO
      IF(LISTBL(N).EQ.NBL) GO TO 523
  524 CONTINUE
      NO = NO + 1
  525 LISTBL(NO)=NBL
      NPL = NFI+1
      NMI = NFI-1
      IF(NFI.EQ. 0) NMI = 83
      IF(NFI.EQ.83) NPL =  0
      LISTBL(NO+1) =NZ+32*NPL
      LISTBL(NO+2) =NZ+32*NMI
      NO = NO + 2
  523 RR(1)=RR(1)+XL*PVX
      RR(2)=RR(2)+XL*PVY
      RR(3)=RR(3)+XL*PVZ
C
C      CHECK WHETHER PARTICLE IS STILL IN THE LEAD GLASS
C
      R1=RR(1)**2+RR(2)**2
      IF(R1.GT.OUTR) GO TO 530
      IF(ZLGPL.LT.RR(3).OR.RR(3).LT.ZLGMI) GO TO 530
C
      GO TO 521
C
C     CALCULATE BLOCK ENERGIES
C
  530 IPART = 0
      VAR(1)= ATAN2(R(2),R(1))
      VAR(2)= R(3)
      IF(LISTBL(1).LT.0) RETURN 1
C ****    EXPEC6 USES BLOCK COUNTING STARTING AT 0
      DO 532 N=1,NO
  532 LISTBL(N) = LISTBL(N) - 1
      CALL EXPEC6(IPART,ENERGY,VAR,DC,NO,LISTBL,BLCFRC,AVGSUM)
C     IF(ESEEN.GT.0.) READCR = ENERGY/ESEEN
C *** CORRECT FOR READOUT THRESHOLD
C
C  AVGSUM NORMS  THE BLCFRC SUM, NEEDED TO KEEP THE RESOLUTION
C  CHANGE 13.4.87, D.PITZL, J.O.
C     ESEEN = 0.
C     DO 533 NB=1,NO
C     EBLOCK = ENERGY*BLCFRC(NB)
C     IF(EBLOCK.GT.28.) ESEEN = ESEEN + EBLOCK
C     ESEEN = ESEEN + EBLOCK
C 533 CONTINUE
      READCR = 1.
C     IF(ESEEN.GT.0.) READCR = ENERGY/ESEEN
C  CHANGED 3.7.86  19.30    J.OLSSON
C     IF(ENERGY.GT.0.) READCR = ESEEN/ENERGY
C ***  UPDATE BLOCK ENERGIES
      DO 531 NB=1,NO
      NBL = LISTBL(NB) + 1
C     AMPLG(NBL) = AMPLG(NBL) + READCR*ENERGY*BLCFRC(NB)
      IF(ABS(AVGSUM).GT..00001)
     $ AMPLG(NBL) = AMPLG(NBL) + READCR*ENERGY*BLCFRC(NB)/AVGSUM
  531 CONTINUE
      R(1) = RR(1)
      R(2) = RR(2)
      R(3) = RR(3)
      RETURN
C
C     END CAP LEAD GLASS
C     FOLLOW PHOTON TO THE FACE OF THE END CAP COUNTERS
C
  700 CONTINUE
      RABS=0.
      DO 701 I3=1,3
  701 RABS=RABS+R(I3)**2
      RABS=SQRT(RABS)
      DO 702 I3=1,3
  702 R(I3)=R(I3)/RABS
      DISTZ=ZENDPL
      IF(R(3).LT.0) DISTZ=ABS(ZENDMI)
      ZLAM=ABS(DISTZ/R(3))
      DO 703 I3=1,3
  703 R(I3)=ZLAM*R(I3)
C
C      SET UP BLOCK ARRAY OF BLOCKS WHICH MIGHT RECIEVE ENERGY
C      FROM THE EL.MG. SHOWER
C
      NO = 0
      RR(3) = R(3)
      DO 710 IC=1,3
      RR(1) = R(1) + (IC-2)*100.
      DO 711 JC=1,3
      RR(2) = R(2) + (JC-2)*100.
      CALL ENDCLG(RR,NBL,*711)
      IF(NO.EQ.0) GO TO 714
      DO 712 N=1,NO
      IF(LISTBL(N).EQ.NBL) GO TO 711
  712 CONTINUE
  714 NO = NO+1
      LISTBL(NO)=NBL
  711 CONTINUE
  710 CONTINUE
C
C     CALCULATE BLOCK ENERGIES
C
  730 IPART = 1
      IF(R(3).LT.0.) IPART = -1
      VAR(1)= R(1)
      VAR(2)= R(2)
      IF(LISTBL(1).LT.0) RETURN 1
C *****     EXPEC6 USES USUAL BLOCK COUNTING STARTING FROM 0
      DO 732 N=1,NO
  732 LISTBL(N) = LISTBL(N) - 1
      CALL EXPEC6(IPART,ENERGY,VAR,DC,NO,LISTBL,BLCFRC,AVGSUM)
C *** CORRECT FOR READOUT THRESHOLD
C     ESEEN = 0.
C     DO 733 NB=1,NO
C     EBLOCK = ENERGY*BLCFRC(NB)
C     IF(EBLOCK.GT.28.) ESEEN = ESEEN + EBLOCK
C 733 CONTINUE
      READCR = 1.
C     IF(ESEEN.GT.0.) READCR = ENERGY/ESEEN
C  CHANGED 3.7.86  19.30    J.OLSSON
C ***  UPDATE BLOCK ENERGIES
      DO 731 NB=1,NO
      NBL = LISTBL(NB) + 1
      IF(ABS(AVGSUM).GT..00001)
     $ AMPLG(NBL) = AMPLG(NBL) + READCR*ENERGY*BLCFRC(NB)/AVGSUM
  731 CONTINUE
C
C UPDATE POSITION VECTOR TO BACK END OF LG, 20 % EXTRA FOR SECURITY
C
      TRKL=DEPEND*1.2/ABS(PVZ)
C
      R(1)=R(1)+PVX*TRKL
      R(2)=R(2)+PVY*TRKL
      R(3)=R(3)+PVZ*TRKL
C
      RETURN
      END
      SUBROUTINE EXPEC6(IPART,E,VAR,DIRC,NBLOCK,IBLIST,BLCFRC,AVGSUM)
C
C     CALCULATE EXPECTED FRACTIONS OF ELECTROMAGNETIC SHOWER
C     IN NBLOCK BLOCKS GIVEN IN A BLOCKLIST IBLIST
C     DISTINGUISH SF5 AND SF6 BLOCKS
C                                   VARIABLES: ...5, ...6
C
C         >>>>>>>>>  MONTE CARLO VERSION  <<<<<<<<<
C
C
C     I/P : IPART        DETECTOR PART -1,0,+1
C           E            SHOWER ENERGY (DEFAULT 500 MEV)
C           VAR          IMPACT POINTS ON LEAD-GLASS SURFACE
C                        VAR(1),VAR(2) = PHI,Z IN BARREL CASE
C                        VAR(1),VAR(2) =   X,Y IN ENDCAP CASE
C           DIRC         DIRECTION COSINE OF SHOWER AXIS (DX,DY,DZ)
C           NBLOCK       LENGTH OF BLOCKLIST
C           IBLIST       LIST OF BLOCKS FOR WHICH THE FRACTIONS
C                        SHOULD BE CALCULATED
C
C     O/P : BLCFRC       FRACTION IN EACH REQUESTED BLOCK *
C                        CORRECTION FOR OUTSIDE FRACTION
C           AVGSUM       ENERGY FRACTION IN BLOCKS CONTAINED IN IBLIST
C
C                                       06/08/82
C               EXTENDED  N.MAGNUSSEN   20/08/85
      EXTERNAL SF5
      EXTERNAL SF6
      EXTERNAL SF65
      EXTERNAL SF56
C
C-----------------------------------------------------------------------
C                            MACRO CGEO1 .... JADE GEOMETRY
C-----------------------------------------------------------------------
C
      COMMON / CGEO1 / BKGAUS,
     +                 RPIP,DRPIP,XRLPIP,   RBPC,DRBPC,XRLBPC,
     +                 RITNK,DRITNK,XRLTKI, R0ROH,DR0ROH,XR0ROH,
     +                 R1ROH,DR1ROH,XR1ROH, R2ROH,DR2ROH,XR2ROH,
     +                 R3ROH,DR3ROH,XR3ROH, ROTNK,DROTNK,XRLTKO,
     +                 RTOF,DRTOF,XRTOF,    RCOIL,DRCOIL,XRCOIL,
     +                 ZJM,DZJM,XRZJM,ZJP,DZJP,XRZJP,ZTKM,DZTKM,XRZTKM,
     +                 ZTKP,DZTKP,XRZTKP,ZBPPL,ZBPMI,ZTOFPL,ZTOFMI,
     +                 XRJETC,RLG,ZLGPL,ZLGMI,OUTR2,CTLIMP,
     +                 CTLIMM,DELFI,BLXY,BLZ,BLDEP,ZENDPL,ZENDMI,DEPEND,
     +                 XHOL1,XHOL2,YHOL1,YHOL2,BLFI
C
C------------------------- END OF MACRO CGEO1 --------------------------
C
C
      COMMON /TRAFO/ STR,CTR,SPR,CPR,CTRCPR,CTRSPR,STRCPR,STRSPR,
     *               JPART,XMOV5,YMOV5,ZMOV5,XMOV6,YMOV6,ZMOV6
      COMMON /PARAME/ P(6,2),MEDNOW
      COMMON / CJTRIG / PI,TWOPI
      DIMENSION VAR(2),DIRC(3),X0(2)
      DIMENSION EL(6),AL5(6),ALL5(6),BL5(6),ED(5),AD5(5),ALD5(5),BD5(5)
      DIMENSION AL6(6),ALL6(6),BL6(6),AD6(5),ALD6(5),BD6(5)
      DIMENSION XL(3),XU(3)
      DIMENSION IBLIST(1),BLCFRC(1)
C
C     SHOWER FUNCTION PARAMETERS ACCORDING TO EGS3 - MONTE CARLO
C
      DATA EL     /100.   ,300.   ,500.   ,700.   ,1000.  ,5000.  /
      DATA AL5    /  3.25 ,  4.59 ,  5.80 ,  6.06 ,  6.72 ,  8.63 /
      DATA ALL5   /   .73 ,  1.45 ,  1.65 ,  1.84 ,  2.03 ,  2.74 /
      DATA BL5    /   .395,   .443,   .447,   .449,   .460,   .455/
      DATA AL6    /  3.25 ,  4.48 ,  5.54 ,  5.87 ,  6.48 ,  8.06 /
      DATA ALL6   /   .73 ,  1.45 ,  1.65 ,  1.84 ,  2.03 ,  2.74 /
      DATA BL6    /   .381,   .429,   .432,   .437,   .447,   .441/
C
      DATA ED     /100.   ,200.   ,200.   ,300.   ,4000.   /
      DATA AD5    /  1.34 ,  1.21 ,  0.26 ,   .66 ,   1.91 /
      DATA ALD5   /   .72 ,   .2  ,   .19 ,   .19 ,    .71 /
      DATA BD5    /   .048,   .004,   .002,   .011,   -.005/
      DATA AD6    /  1.23 ,  1.06 ,   .33 ,   .61 ,   1.58 /
      DATA ALD6   /   .72 ,   .2  ,   .19 ,   .19 ,    .71 /
      DATA BD6    /   .048,   .003,   .005,   .010,   -.006/
C
      DATA FACMC /.9013/
      DATA X0/23.90,16.60/
C     DATA ICO/0/
      DATA ZSF6/318./
C
      JPART = IPART
      DO 116 IBL = 1,NBLOCK
      BLCFRC(IBL) = 0.
  116 CONTINUE
      AVGSUM = 0.
C
C     IMPACT POINT ON LG SURFACE ??
C
      IF(ICHECK(VAR(1),VAR(2),IPART).NE.1) RETURN
C
C     TRANSFORMATION PARAMETERS
C
      DX = DIRC(1)
      DY = DIRC(2)
      DZ = DIRC(3)
      IF(IPART.EQ.0) GOTO 1
      ZENDCP = IPART*ZENDPL
      XMOV = VAR(1)
      YMOV = VAR(2)
      ZMOV = ZENDCP
      GOTO 2
    1 XMOV = RLG*COS(VAR(1))
      YMOV = RLG*SIN(VAR(1))
      ZMOV = VAR(2)
    2 CTR = DZ/SQRT(DX*DX+DY*DY+DZ*DZ)
      STR = SQRT(1.-CTR*CTR)
      PR=ATAN2(DY,DX)
      IF(PR.LT.0.) PR=PR+TWOPI
      SPR = SIN(PR)
      CPR = COS(PR)
      CTRCPR = CTR*CPR
      CTRSPR = CTR*SPR
      STRCPR = STR*CPR
      STRSPR = STR*SPR
C
C     CHOOSE SHOWER START POINT (THIS IS THE MC-VERSION !)
C
      SSPOI5 = -X0(1)*ALOG(RN(DUMMY))/.773278
      SSPOI6 = X0(2)*SSPOI5/X0(1)
      XMOV5 = XMOV + SSPOI5*DX
      YMOV5 = YMOV + SSPOI5*DY
      ZMOV5 = ZMOV + SSPOI5*DZ
      XMOV6 = XMOV + SSPOI6*DX
      YMOV6 = YMOV + SSPOI6*DY
      ZMOV6 = ZMOV + SSPOI6*DZ
C
C     SHOWER PARAMETERS (PLO'I'5,6  INTERPOLATED FROM PARAM.)
C
      IF(E.LE.0.) E=500.
      IF(E.LT.5000.) GOTO 3
      PLO15 = AL5(6)
      PLO25 = ALL5(6)
      PLO35 = BL5(6)
      PLO16 = AL6(6)
      PLO26 = ALL6(6)
      PLO36 = BL6(6)
      GOTO 9
    3 IF(E.LT.1000.) GOTO 4
      EFAC = (E-EL(5))/ED(5)
      PLO15 = AL5(5)+AD5(5)*EFAC
      PLO25 = ALL5(5)+ALD5(5)*EFAC
      PLO35 = BL5(5)+BD5(5)*EFAC
      PLO16 = AL6(5)+AD6(5)*EFAC
      PLO26 = ALL6(5)+ALD6(5)*EFAC
      PLO36 = BL6(5)+BD6(5)*EFAC
      GOTO 9
    4 IF(E.LT.700.)  GOTO 5
      EFAC = (E-EL(4))/ED(4)
      PLO15 = AL5(4)+AD5(4)*EFAC
      PLO25 = ALL5(4)+ALD5(4)*EFAC
      PLO35 = BL5(4)+BD5(4)*EFAC
      PLO16 = AL6(4)+AD6(4)*EFAC
      PLO26 = ALL6(4)+ALD6(4)*EFAC
      PLO36 = BL6(4)+BD6(4)*EFAC
      GOTO 9
    5 IF(E.LT.500.)  GOTO 6
      EFAC = (E-EL(3))/ED(3)
      PLO15 = AL5(3)+AD5(3)*EFAC
      PLO25 = ALL5(3)+ALD5(3)*EFAC
      PLO35 = BL5(3)+BD5(3)*EFAC
      PLO16 = AL6(3)+AD6(3)*EFAC
      PLO26 = ALL6(3)+ALD6(3)*EFAC
      PLO36 = BL6(3)+BD6(3)*EFAC
      GOTO 9
    6 IF(E.LT.300.)  GOTO 7
      EFAC = (E-EL(2))/ED(2)
      PLO15 = AL5(2)+AD5(2)*EFAC
      PLO25 = ALL5(2)+ALD5(2)*EFAC
      PLO35 = BL5(2)+BD5(2)*EFAC
      PLO16 = AL6(2)+AD6(2)*EFAC
      PLO26 = ALL6(2)+ALD6(2)*EFAC
      PLO36 = BL6(2)+BD6(2)*EFAC
      GOTO 9
    7 IF(E.LT.100.)  GOTO 8
      EFAC = (E-EL(1))/ED(1)
      PLO15 = AL5(1)+AD5(1)*EFAC
      PLO25 = ALL5(1)+ALD5(1)*EFAC
      PLO35 = BL5(1)+BD5(1)*EFAC
      PLO16 = AL6(1)+AD6(1)*EFAC
      PLO26 = ALL6(1)+ALD6(1)*EFAC
      PLO36 = BL6(1)+BD6(1)*EFAC
      GOTO 9
    8 PLO15 = AL5(1)
      PLO25 = ALL5(1)
      PLO35 = BL5(1)
      PLO16 = AL6(1)
      PLO26 = ALL6(1)
      PLO36 = BL6(1)
C
    9 P(1,1) = PLO25
      P(2,1) = PLO15*3.293
      P(3,1) = PLO35
      P(4,1) = 3.134
      P(5,1) = PLO15*.267
      P(6,1) = .6354
      P(1,2) = PLO26
      P(2,2) = PLO16*2.832
      P(3,2) = PLO36
      P(4,2) = 3.017
      P(5,2) = PLO16*.266
      P(6,2) = .60333
C
C     NORMALIZE INTEGRALS
C
      P2A15=PLO25+1.
      P2A16=PLO26+1.
C     FNORM5 = FACMC*60500.0*PLO15*GAMMA(P2A15)/PLO35**P2A15
C     FNORM6 = FACMC*22000.0*PLO16*GAMMA(P2A16)/PLO36**P2A16
C     FNORM5 = FACMC*59800.0*PLO15*GAMMA(P2A15)/PLO35**P2A15
C     FNORM6 = FACMC*21800.0*PLO16*GAMMA(P2A16)/PLO36**P2A16
C     FNORM5 = FACMC*62618.0*PLO15*GAMMA(P2A15)/PLO35**P2A15
C     FNORM6 = FACMC*21936.0*PLO16*GAMMA(P2A16)/PLO36**P2A16
      FNORM5 = FACMC*59000.0*PLO15*GAMMA(P2A15)/PLO35**P2A15
      FNORM6 = FACMC*21220.0*PLO16*GAMMA(P2A16)/PLO36**P2A16
C
C     INTEGRATION
C
C
C     CHECK TOTAL VALUE ON REQUEST
C
      IF(NBLOCK.GT.0) GOTO 113
      XL(1) = RLG
      IF(IPART.NE.0) XL(1) = VAR(1)-500.
      XL(2) = VAR(1)-2.
      IF(IPART.NE.0) XL(2) = VAR(2)-500.
      XL(3) = VAR(2)-500.
      IF(IPART.EQ. 1) XL(3) = ZENDPL
      IF(IPART.EQ.-1) XL(3) = ZENDMI-1000.
      XU(1) = RLG+1000.
      IF(IPART.NE.0) XU(1) = VAR(1)+500.
      XU(2) = VAR(1)+2.
      IF(IPART.NE.0) XU(2) = VAR(2)+500.
      XU(3) = VAR(2)+500.
      IF(IPART.EQ. 1) XU(3) = ZENDPL+1000.
      IF(IPART.EQ.-1) XU(3) = ZENDMI
      CALL WEGAS(SF5,.001,1000,9,XL,XU,AVG,SD,C2)
      AVG = AVG/FNORM5
      SD = SD/FNORM5
C     WRITE(6,1133) AVG,SD,C2
C1133 FORMAT(10X,'---> WEGAS INTEGRAL RESULT : ',F12.5,/,
C    *       10X,'---> ERROR + CHI**2 VALUE  : ',2F12.5)
      RETURN
C
C     PERFORM INTEGRATION OVER REQUESTED BLOCKS
C
  113 MED5=0
      MED6=0
      BEHIN5=E*X0(1)/100.
      BEHIN6=E*X0(2)/100.
      BFRONT=10.
      DO 1131 IB=1,NBLOCK
       IBLOCK = IBLIST(IB)
       IWHICH = 5
       J=MOD(IBLOCK,32)
       IF(J.GE.13.AND.J.LE.18.AND.IPART.EQ.0) IWHICH = 6
       MED = IWHICH
       IF(MED.EQ.5) MED5=1
       IF(MED.EQ.6) MED6=1
 1131 CONTINUE
      IF(MED5.EQ.1.AND.MED6.EQ.0) GOTO 1145
      IF(MED5.EQ.0.AND.MED6.EQ.1) GOTO 1146
      IF(MED5.EQ.1.AND.MED6.EQ.1) GOTO 1147
C
C  ONLY SF5-BLOCKS IN CLUSTER
C
1145  DO 1141 IBL = 1,NBLOCK
       IBLOCK = IBLIST(IBL)
       CALL LIMINT(IBLOCK,IPART,IER,XL(1),XL(2),XL(3),XU(1),XU(2),XU(3))
       IF(IER.GT.0) GOTO 1141
       IF(IPART.EQ. 0) XU(1) = XU(1) + BEHIN5
       IF(IPART.EQ. 1) XU(3) = XU(3) + BEHIN5
       IF(IPART.EQ.-1) XL(3) = XL(3) - BEHIN5
       IF(IPART.EQ. 0) XL(1) = XL(1) - BFRONT
       IF(IPART.EQ. 1) XL(3) = XL(3) - BFRONT
       IF(IPART.EQ.-1) XU(3) = XU(3) + BFRONT
C
C     INTEGRATION OVER A SF5-BLOCK
C
       CALL WEGAS(SF5,.01,200,3,XL,XU,AVG,SD,C2)
       AVG = AVG/FNORM5
       SD = SD/FNORM5
       BLCFRC(IBL) = AVG
       AVGSUM = AVGSUM+AVG
1141  CONTINUE
      GOTO 115
C
C   ONLY SF6-BLOCKS IN CLUSTER
C
1146  CONTINUE
      DO 1142 IBL=1,NBLOCK
       IBLOCK = IBLIST(IBL)
       CALL LIMINT(IBLOCK,IPART,IER,XL(1),XL(2),XL(3),XU(1),XU(2),XU(3))
       IF(IER.GT.0) GOTO 1142
       XU(1) = XU(1) + BEHIN6
       XL(1) = XL(1) - BFRONT
       CALL WEGAS(SF6,.01,200,3,XL,XU,AVG,SD,C2)
       AVG = AVG/FNORM6
       SD = SD/FNORM6
C     IF (ICO.GT.1000) GOTO 2135
C     ICO=ICO+1
C     ESEE=E*AVG
C     WRITE(6,2134) IBLOCK,AVG,SD,C2,ESEE,XMOV6,YMOV6,ZMOV6
C2134 FORMAT(10X,'*--> WEGAS INTEGRAL RESULT SF6#: ',I8,F12.5,/,
C    *       10X,'*--> ERROR + CHI**2 VALUE+ESEEN : ',3F12.5,/,
C    *       10X,'*--> X,Y,ZMOV6: ',3F12.5)
C2135 CONTINUE
       BLCFRC(IBL) = AVG
       AVGSUM = AVGSUM+AVG
1142  CONTINUE
      GOTO 115
C
C     BORDER CROSSING; SF5*SF6 IN CLUSTER
C
1147  CONTINUE
      MEDSRT=5
      IF(ABS(ZMOV6).LT.ZSF6) MEDSRT=6
      DO 1143 IBL=1,NBLOCK
       IBLOCK = IBLIST(IBL)
       CALL LIMINT(IBLOCK,IPART,IER,XL(1),XL(2),XL(3),XU(1),XU(2),XU(3))
       IF(IER.GT.0) GOTO 1143
       IWHICH = 5
       JJ=MOD(IBLOCK,32)
       IF(JJ.GE.13.AND.JJ.LE.18.AND.IPART.EQ.0) IWHICH = 6
       MEDNOW=IWHICH
       IF(MEDSRT.EQ.6.AND.MEDNOW.EQ.5) GOTO 1159
       IF(MEDSRT.EQ.5.AND.MEDNOW.EQ.6) GOTO 1169
       IF(MEDSRT.EQ.5.AND.MEDNOW.EQ.5) GOTO 1149
C START IN SF6 AND NOW IN SF6
       XU(1) = XU(1) + BEHIN6
       XL(1) = XL(1) - BFRONT
        CALL WEGAS(SF6,.01,200,3,XL,XU,AVG,SD,C2)
       AVG = AVG/FNORM6
       SD = SD/FNORM6
C     IF (ICO.GT.1000) GOTO 1135
C     ICO=ICO+1
C     ESEE=E*AVG
C     WRITE(6,1134) IBLOCK,AVG,SD,C2,ESEE
C1134 FORMAT(10X,'---> WEGAS INTEGRAL RESULT SF6#: ',I8,F12.5,/,
C    *       10X,'---> ERROR + CHI**2 VALUE+ESEEN : ',3F12.5)
C1135 CONTINUE
       BLCFRC(IBL) = AVG
       AVGSUM = AVGSUM+AVG
       GOTO 1143
C START IN SF5 AND NOW IN SF5
1149   XU(1) = XU(1) + BEHIN5
       XL(1) = XL(1) - BFRONT
        CALL WEGAS(SF5,.01,200,3,XL,XU,AVG,SD,C2)
       AVG = AVG/FNORM5
       SD = SD/FNORM5
C     IF (ICO.GT.1000) GOTO 1335
C     ICO=ICO+1
C     ESEE=E*AVG
C     WRITE(6,1334) IBLOCK,AVG,SD,C2,ESEE
C1334 FORMAT(10X,'---> WEGAS INTEGRAL RESULT SF5#: ',I8,F12.5,/,
C    *       10X,'---> ERROR + CHI**2 VALUE+ESEEN : ',3F12.5)
C1335 CONTINUE
       BLCFRC(IBL) = AVG
       AVGSUM = AVGSUM+AVG
       GOTO 1143
C START IN SF6 AND NOW IN SF5
1159  XU(1) = XU(1) + BEHIN5
      XL(1) = XL(1) - BFRONT
      CALL WEGAS(SF65,.01,200,3,XL,XU,AVG,SD,C2)
      AVG = AVG/FNORM5
      SD = SD/FNORM5
C     IF (ICO.GT.1000) GOTO 1137
C     ICO=ICO+1
C     ESEE=E*AVG
C     WRITE(6,1136) IBLOCK,AVG,SD,C2,ESEE
C1136 FORMAT(10X,'---> WEGAS INTEGRAL RESULT SF65#: ',I8,F12.5,/,
C    *       10X,'---> ERROR + CHI**2 VALUE+ESEE  : ',3F12.5)
C1137 CONTINUE
      BLCFRC(IBL) = AVG
      AVGSUM = AVGSUM+AVG
      GOTO 1143
C START IN SF5 AND NOW IN SF6
1169  XU(1) = XU(1) + BEHIN6
      XL(1) = XL(1) - BFRONT
      CALL WEGAS(SF56,.01,200,3,XL,XU,AVG,SD,C2)
      AVG = AVG/FNORM6
      SD = SD/FNORM6
C     IF (ICO.GT.1000) GOTO 1137
C     ICO=ICO+1
C     ESEE=E*AVG
C     WRITE(6,1336) IBLOCK,AVG,SD,C2,ESEE
C1336 FORMAT(10X,'---> WEGAS INTEGRAL RESULT SF56#: ',I8,F12.5,/,
C    *       10X,'---> ERROR + CHI**2 VALUE+ESEE  : ',3F12.5)
C1137 CONTINUE
      BLCFRC(IBL) = AVG
      AVGSUM = AVGSUM+AVG
1143  CONTINUE
 115  CONTINUE
C
      AVGDIF = 1.-AVGSUM
C
      RETURN
      END
C
C
      REAL FUNCTION SF5(X)
C
C     FUNCTION TO BE INTEGRATED BY WEGAS
C                                       16/06/82
C                     CHANGED   N.M.    26/08/85
C        CHANGED TO MC VERSION  N.M.    11/10/85
C
      COMMON /TRAFO/ STR,CTR,SPR,CPR,CTRCPR,CTRSPR,STRCPR,STRSPR,
     *               JPART,XMOV5,YMOV5,ZMOV5,XMOV6,YMOV6,ZMOV6
      COMMON /PARAME/ P(6,2),MEDNOW
      REAL X(3)
      DATA X0 /23.90/
C
      SF5 = 0.
C
C     TRANSFORM INTO CARTESIAN DETECTOR SYSTEM
C
      IF(JPART.NE.0) GOTO 1
      XDET = X(1)*COS(X(2))
      YDET = X(1)*SIN(X(2))
      ZDET = X(3)
      FJACOB = X(1)
      GOTO 2
    1 XDET = X(1)
      YDET = X(2)
      ZDET = X(3)
      FJACOB = 1.
C
C     LINEAR TRANSFORMATION TO SHOWER START POINT
C
    2 XLTR = XDET-XMOV5
      YLTR = YDET-YMOV5
      ZLTR = ZDET-ZMOV5
C
C     ROTATE INTO SHOWER SYSTEM (FIRST CHECK SIGN OF ZSHW)
C
      ZSHW  =  XLTR * STRCPR   +   YLTR * STRSPR   +   ZLTR * CTR
C
      IF(ZSHW.LE.0.) RETURN
C
      XSHW  =  XLTR * CTRCPR   +   YLTR * CTRSPR   -   ZLTR * STR
      YSHW  = -XLTR * SPR      +   YLTR * CPR
C
C     SCALE IN RADIATION LENGTHS
C
      RSHS = SQRT(XSHW*XSHW+YSHW*YSHW)/X0
      ZSHS = ZSHW/X0
C
C     COMPRESS SHOWER FUNCTION
C
      P3Z = P(3,1)*ZSHS
      REX = -1.*P(1,1)*ALOG(ZSHS)
      EX1 = -1.*(REX+P3Z+P(4,1)*RSHS)
      EX2 = -1.*(REX+P3Z+P(6,1)*RSHS)
      IF(EX1.LT.-10..OR.EX2.LT.-10.) RETURN
C
      SF5 =  P(2,1)*EXP(EX1) + P(5,1)*EXP(EX2)
C
C     JACOBIAN
C
      SF5 = SF5*FJACOB
C
      RETURN
      END
C
C
      REAL FUNCTION SF6(X)
C
C     FUNCTION TO BE INTEGRATED BY WEGAS
C                                       16/06/82
C                     CHANGED   N.M.    26/08/85
C        CHANGED TO MC VERSION  N.M.    11/10/85
C
      COMMON /TRAFO/ STR,CTR,SPR,CPR,CTRCPR,CTRSPR,STRCPR,STRSPR,
     *               JPART,XMOV5,YMOV5,ZMOV5,XMOV6,YMOV6,ZMOV6
      COMMON /PARAME/ P(6,2),MEDNOW
      REAL X(3)
      DATA X0 /16.60/
C
      SF6 = 0.
C
C     TRANSFORM INTO CARTESIAN DETECTOR SYSTEM
C
      IF(JPART.NE.0) GOTO 1
      XDET = X(1)*COS(X(2))
      YDET = X(1)*SIN(X(2))
      ZDET = X(3)
      FJACOB = X(1)
      GOTO 2
    1 XDET = X(1)
      YDET = X(2)
      ZDET = X(3)
      FJACOB = 1.
C
C     LINEAR TRANSFORMATION TO SHOWER START POINT
C
    2 XLTR = XDET-XMOV6
      YLTR = YDET-YMOV6
      ZLTR = ZDET-ZMOV6
C
C     ROTATE INTO SHOWER SYSTEM (FIRST CHECK SIGN OF ZSHW)
C
      ZSHW  =  XLTR * STRCPR   +   YLTR * STRSPR   +   ZLTR * CTR
C
      IF(ZSHW.LE.0.) RETURN
C
      XSHW  =  XLTR * CTRCPR   +   YLTR * CTRSPR   -   ZLTR * STR
      YSHW  = -XLTR * SPR      +   YLTR * CPR
C
C     SCALE IN RADIATION LENGTHS
C
      RSHS = SQRT(XSHW*XSHW+YSHW*YSHW)/X0
      ZSHS = ZSHW/X0
C
C     COMPRESS SHOWER FUNCTION
C
      P3Z = P(3,2)*ZSHS
      REX = -1.*P(1,2)*ALOG(ZSHS)
      EX1 = -1.*(REX+P3Z+P(4,2)*RSHS)
      EX2 = -1.*(REX+P3Z+P(6,2)*RSHS)
      IF(EX1.LT.-10..OR.EX2.LT.-10.) RETURN
C
      SF6 =  P(2,2)*EXP(EX1) + P(5,2)*EXP(EX2)
C
C     JACOBIAN
C
      SF6 = SF6*FJACOB
C
      RETURN
      END
C
      REAL FUNCTION SF65(X)
C
C     FUNCTION TO BE INTEGRATED BY WEGAS
C     SF6/SF5 BORDER CROSSED FROM SF6 INTO SF5
C                                       16/06/82
C                     CHANGED   N.M.    26/08/85
C        CHANGED TO MC VERSION  N.M.    11/10/85
C
      COMMON /TRAFO/ STR,CTR,SPR,CPR,CTRCPR,CTRSPR,STRCPR,STRSPR,
     *               JPART,XMOV5,YMOV5,ZMOV5,XMOV6,YMOV6,ZMOV6
      COMMON /PARAME/ P(6,2),MEDNOW
      REAL X(3),X0(2)
      DATA X0 /23.90,16.60/,ZSF6/318./
      DATA IC /0/, ICAL65/0/
C
      SF65 = 0.
C
C     TRANSFORM INTO CARTESIAN DETECTOR SYSTEM
C
      IF(JPART.NE.0) GOTO 1
      XDET = X(1)*COS(X(2))
      YDET = X(1)*SIN(X(2))
      ZDET = X(3)
      FJACOB = X(1)
      GOTO 2
    1 XDET = X(1)
      YDET = X(2)
      ZDET = X(3)
      FJACOB = 1.
C
C     LINEAR TRANSFORMATION TO SHOWER START POINT
C     IMPLICTLY ABS(.DET) GT ABS(.MOV)
C
    2 XLTR = XDET-XMOV6
      YLTR = YDET-YMOV6
C     IF(ZMOV6.LT.0.)
C    *WRITE (6,9001) JPART,XDET,YDET,ZDET,XMOV6,YMOV6,ZMOV6
C
C     DETECTOR-COORDINATES OF PUNCHTHROUGH
C     INTEGRATION NOW IN SF5-BLOCK: POSITIVE Z-SIDE
C
      ZLTR = ZDET-ZMOV6
      IF(ZDET.LT.0.) GOTO 200
      DZ1 = ZSF6-ZMOV6
      DZ2 = ZDET-ZSF6
C     ZSHWP = DZ1/CTR
      GOTO 300
C     INTEGRATION NOW IN SF5-BLOCK: NEGATIVE Z-SIDE
 200  DZ1 = ZSF6+ZMOV6
      DZ2 = -ZSF6-ZDET
 300  CONTINUE
      ZSHWP = ABS(DZ1/CTR)
C     IF(ZMOV6.LT.0.)
C    *WRITE (6,9002) XLTR,YLTR,ZLTR,ZSHWP,DZ1,DZ2
C     IF(ZMOV6.LT.0.)
C    *WRITE (6,9004) STR,CTR,SPR,CPR
C
C     ROTATE INTO SHOWER SYSTEM (FIRST CHECK SIGN OF ZSHW)
C
      ZSHW  = XLTR * STRCPR   +  YLTR * STRSPR   +  ZLTR * CTR
C     IF(ZMOV6.LT.0.)
C    *WRITE(6,9005) ZSHW
C
      IF(ZSHW.LE.0.) RETURN
C
      XSHW  =  XLTR * CTRCPR   +   YLTR * CTRSPR   -   ZLTR * STR
      YSHW  = -XLTR * SPR      +   YLTR * CPR
C
C     SCALE IN RADIATION LENGTHS
C
      IF(ZSHWP.LT.0.) ICAL65 = ICAL65 + 1
      IF(ZSHWP.LT.0..AND.ICAL65.LT.100)
     *WRITE (6,9002) XLTR,YLTR,ZLTR,ZSHWP,DZ1,DZ2
      IF(ZSHWP.LT.0..AND.ICAL65.LT.100)
     *WRITE (6,9004) STR,CTR,SPR,CPR
      ZSHWP = ABS(ZSHWP)/X0(2)
      ZSHS = ZSHW/X0(2)
      DELZ = ZSHW - ZSHWP
      IF(DELZ.LE.0.) GOTO 20
      ZSHS = ZSHWP/X0(2) + DELZ/X0(1)
C  DEVIDE RSHW ACCORDING TO DZ1,DZ2
 20   RSHW = SQRT(XSHW*XSHW+YSHW*YSHW)
      RSHW1 = 0.
      IF(ZSHW.GT.ZSHWP) GOTO 25
      RSHW0 = DZ1/STR
      RSHW1 = -ZSHW*CTR/STR
      IF(CTR.LT.0.) RSHW1 = ZSHW*CTR/STR
      RSHW1 = RSHW1 + RSHW0
25    RSHW2 = RSHW - RSHW1
      RSHS = RSHW1/X0(2)+RSHW2/X0(1)
C
C     IF(DELZ.GT.0.) RSHS = RSHS*X0(2)/X0(1)
C     IF(ZMOV6.LT.0.)
      IF(RSHW1.LT.0..OR.RSHW2.LT.0.)
     *WRITE (6,9004) STR,CTR,SPR,CPR
      IF(RSHW1.LT.0..OR.RSHW2.LT.0.)
     *WRITE (6,9002) XLTR,YLTR,ZLTR,ZSHWP,DZ1,DZ2
      IF(RSHW1.LT.0..OR.RSHW2.LT.0.)
     *WRITE(6,9006) XSHW,YSHW,RSHS,ZSHS,RSHW1,RSHW2,RSHW
      IF(MEDNOW.EQ.5) GOTO 999
      IF(IC.GE.100) GOTO 999
C      IF(MEDNOW.NE.5) WRITE (6,6000) MEDNOW
      IC=IC+1
  999 CONTINUE
C
C     COMPRESS SHOWER FUNCTION
C
      P3Z = P(3,1)*ZSHS
      REX = -1.*P(1,1)*ALOG(ZSHS)
      EX1 = -1.*(REX+P3Z+P(4,1)*RSHS)
      EX2 = -1.*(REX+P3Z+P(6,1)*RSHS)
C     IF(ZMOV6.LT.0.)
C    *WRITE(6,9007) EX1,EX2,P3Z,REX,P(3,1),P(1,1)
      IF(EX1.LT.-10..OR.EX2.LT.-10.) RETURN
C
      SF65 =  P(2,1)*EXP(EX1) + P(5,1)*EXP(EX2)
C     IF(ZMOV6.LT.0.)
C    *WRITE(6,9008) SF65

C
C     JACOBIAN
C
 4    SF65 = SF65*FJACOB
 6000 FORMAT(' ***SF65**** MEDNOW UNUSUAL=',I8)
 9001 FORMAT(' JPART,XDET,YDET,ZDET,XMOV6,YMOV6,ZMOV6  ',I6,6F10.4)
 9002 FORMAT(' XLTR,YLTR,ZLTR,ZSHWP,DZ1,DZ2           ',6F12.4)
 9004 FORMAT(' STR,CTR,SPR,CPR      ',4F10.5)
 9005 FORMAT(' ZSHW   ',1F12.5)
 9006 FORMAT(' XSHW,YSHW,RSHS,ZSHS,RSW1,RSHW2,RSHW '/7F10.5)
 9007 FORMAT(' EX1,EX2,P3Z,REX,P(3,2),P(1,2) ',6F12.7)
 9008 FORMAT(' S F 6 5 ',1E12.6)
C
      RETURN
      END
C
      REAL FUNCTION SF56(X)
C
C     FUNCTION TO BE INTEGRATED BY WEGAS
C     SF6/SF5 BORDER CROSSED FROM SF5 INTO SF6
C                                       16/06/82
C                     CHANGED   N.M.    26/08/85
C        CHANGED TO MC VERSION  N.M.    11/10/85
C        CHANGED AGAIN 2.2.87          N.M. / J.O.
C
      COMMON /TRAFO/ STR,CTR,SPR,CPR,CTRCPR,CTRSPR,STRCPR,STRSPR,
     *               JPART,XMOV5,YMOV5,ZMOV5,XMOV6,YMOV6,ZMOV6
      COMMON /PARAME/ P(6,2),MEDNOW
      REAL X(3),X0(2)
      DATA X0 /23.90,16.60/,ZSF6/318./
      DATA IC /0/, ICAL56/0/
C
      SF56 = 0.
C
C     TRANSFORM INTO CARTESIAN DETECTOR SYSTEM
C
      IF(JPART.NE.0) GOTO 1
      XDET = X(1)*COS(X(2))
      YDET = X(1)*SIN(X(2))
      ZDET = X(3)
      FJACOB = X(1)
      GOTO 2
    1 XDET = X(1)
      YDET = X(2)
      ZDET = X(3)
      FJACOB = 1.
C
C     LINEAR TRANSFORMATION TO SHOWER START POINT
C     IMPLICTLY ABS(.DET) GT ABS(.MOV)
C
    2 XLTR = XDET-XMOV5
      YLTR = YDET-YMOV5
C     WRITE (6,9001) JPART,XDET,YDET,ZDET,XMOV5,YMOV5,ZMOV5
C
C     DETECTOR-COORDINATES OF PUNCHTHROUGH
C     INTEGRATION NOW IN SF5-BLOCK: POSITIVE Z-SIDE
C
      ZLTR = ZDET-ZMOV5
      IF(ZDET.LT.0.) GOTO 200
      DZ1 = ZMOV5-ZSF6
      DZ2 = ZSF6 - ZDET
      GOTO 300
 200  DZ1 = -ZSF6-ZMOV5
      DZ2 = ZDET + ZSF6
 300  CONTINUE
C     WRITE (6,9002) XLTR,YLTR,ZLTR
C     WRITE (6,9004) STR,CTR,SPR,CPR
C
C     ROTATE INTO SHOWER SYSTEM (FIRST CHECK SIGN OF ZSHW)
C
      ZSHW  = XLTR * STRCPR   +  YLTR * STRSPR   +  ZLTR * CTR
C     WRITE(6,9005) ZSHW
C
      IF(ZSHW.LE.0.) RETURN
C
      XSHW  =  XLTR * CTRCPR   +   YLTR * CTRSPR   -   ZLTR * STR
      YSHW  = -XLTR * SPR      +   YLTR * CPR
C
C     SCALE IN RADIATION LENGTHS
C
      ZSHS = ZSHW/X0(1)
C  DEVIDE RSHW ACCORDING TO DZ1,DZ2
      RSHW = SQRT(XSHW*XSHW+YSHW*YSHW)
      RSHW1 = RSHW*ABS(DZ1)/(ABS(DZ1)+ABS(DZ2))
      RSHW2 = RSHW - RSHW1
      RSHS = RSHW1/X0(1)+RSHW2/X0(2)
C
      IF(RSHW1.LT.0..OR.RSHW2.LT.0.) ICAL56 = ICAL56 + 1
      IF(ICAL56.GT.100) GO TO 1304
C
      IF(RSHW1.LT.0..OR.RSHW2.LT.0.)
     *WRITE (6,9001) JPART,XDET,YDET,ZDET,XMOV5,YMOV5,ZMOV5
      IF(RSHW1.LT.0..OR.RSHW2.LT.0.)
     *WRITE (6,9004) STR,CTR,SPR,CPR
      IF(RSHW1.LT.0..OR.RSHW2.LT.0.)
     *WRITE (6,9002) XLTR,YLTR,ZLTR,ZSHWP,DZ1,DZ2
      IF(RSHW1.LT.0..OR.RSHW2.LT.0.)
     *WRITE(6,9006) XSHW,YSHW,RSHS,ZSHS,RSHW1,RSHW2,RSHW
C
1304  CONTINUE
      IF(MEDNOW.EQ.6) GOTO 999
      IF(IC.GE.100) GOTO 999
       IF(MEDNOW.NE.6) WRITE (6,6000) MEDNOW
      IC=IC+1
  999 CONTINUE
C
C     COMPRESS SHOWER FUNCTION
C
      P3Z = P(3,2)*ZSHS
      REX = -1.*P(1,2)*ALOG(ZSHS)
      EX1 = -1.*(REX+P3Z+P(4,2)*RSHS)
      EX2 = -1.*(REX+P3Z+P(6,2)*RSHS)
C     WRITE(6,9007) EX1,EX2,P3Z,REX,P(3,2),P(1,2)
      IF(EX1.LT.-10..OR.EX2.LT.-10.) RETURN
C
      SF56 =  P(2,2)*EXP(EX1) + P(5,2)*EXP(EX2)
C     WRITE(6,9008) SF56

C
C     JACOBIAN
C
 4    SF56 = SF56*FJACOB
 9001 FORMAT(' JPART,XDET,YDET,ZDET,XMOV6,YMOV6,ZMOV6  ',I6,6F10.4)
 9002 FORMAT(' XLTR,YLTR,ZLTR            ',3F12.4)
 9004 FORMAT(' STR,CTR,SPR,CPR      ',4F10.5)
 9005 FORMAT(' ZSHW   ',1F12.5)
 9006 FORMAT(' XSHW,YSHW,RSHS,ZSHS,RSHW1,RSHW2,RSHW '/7F10.5)
 9007 FORMAT(' EX1,EX2,P3Z,REX,P(3,2),P(1,2) ',6F12.7)
 9008 FORMAT(' S F 5 6 ',1E12.6)
 6000 FORMAT(' ***SF56**** MEDNOW UNUSUAL=',I8)
C
      RETURN
      END
