CPV/cpvrec.F

   1       SUBROUTINE CPVREC(DzCPV,DyCPV,DxCPV,Nhit,Xin,Pin,Nout,Xout)
   2 C-----------------------------------------------------------------------
   3 C     Reconstruction routing for the CPV detector
   4 C
   5 C     Author: Serguei Sadovski, IHEP, Protvino
   6 C     Created:       25 December 1998
   7 C     Last modified: 17 September 1999
   8 C-----------------------------------------------------------------------
   9 C
  10 C CPV Frame:  Z coordinate - along beam direction
  11 C ==========  Y coordinate - transverse beam direction in the CPV plane
  12 C             X coordinate - transverse CPV plane, positive direction - along
  13 C                            particles going from the interaction point
  14 C            (0,0) in the center of the CPV plane
  15 C
  16 C Input Parameters:
  17 C =================
  18 C DzCPV     - half width (cm) of CPV plane along beam
  19 C DyCPV     - half width (cm) of CPV plane transverse beam
  20 C DxCPV     - half width (cm) of CPV
  21 C Nhit      - number of particles in the front of CPV
  22 C
  23 C Xin(1,i)  - Z coordinate (cm) of i-hit in CPV frame (0,0 in the CPV Center)
  24 C Xin(2,i)  - Y coordinate (cm) of i-hit in CPV frame (0,0 in the CPV Center)
  25 C
  26 C Pin(1,i)  - Pz
  27 C Pin(2,i)  - Py      for particle i
  28 C Pin(3,i)  - Px
  29 C Pin(4,i)  - dE/dX
  30 C
  31 C Output Parameters - reconstructed coordinates of particles in the CPV:
  32 C =================
  33 C Nout      - number of the reconstructed hits (<5000) in CPV plane,
  34 C             if (Nout.le.0) - error output: reconstruction fails ...
  35 C
  36 C Xout(1,j) - Z coordinate (cm) of j-hit in CPV frame (0,0 in the CPV Center)
  37 C Xout(2,j) - Y coordinate (cm) of j-hit in CPV frame (0,0 in the CPV Center)
  38 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  39       DIMENSION Xin(2,Nhit),Pin(4,Nhit),Xout(2,5000)             ! Inp/Out
  40 C
  41       COMMON/GAMMA/ NGAM,EGAM(5000),XGAM(5000),YGAM(5000),CHGAM(5000)! From Rec
  42       COMMON/CONSTA/ NcelZ,NcelY,IDELT,IDELA,CelZ,CelY
  43 C
  44       REAL   *4 EGM,ZGM,YGM, ZYEgm(3,5)
  45       REAL   *8 Pp
  46       INTEGER*2               ICPV(104,176)                      ! CPV 104x176
  47       DATA   NcelZ,NcelY,CelZ,CelY/104,176, 2.26, 1.13 /         ! in cm
  48       DATA   NgamZ,NgamY,Chit,Cnoise/5,9  , 1000., 30. /
  49       DATA   CelWr,RhoWr,IDELT,IDELA/ 0.5   , .01, 90,  90 /     ! 09.01.99
  50       DATA   Detr / 0.1 /  !  Relative energy fluctuation in Track for 100 e-
  51 C
  52       Nout = 0
  53       Ncpv = 0
  54       IF (Nhit.eq.0) RETURN
  55 C!-                              - - -  C O N S T A N T S  - - -
  56       ECPV = 0.
  57       NZMX = 0
  58       NZMN = NcelZ
  59       NYMX = 0
  60       NYMN = NcelY
  61 C
  62       CALL VZERO(ICPV,NcelZ*NcelY/2)
  63 C
  64       DO 1000 Ihit=1,Nhit
  65 C
  66       Pp   =  SQRT( Pin(1,Ihit)**2 +  Pin(2,Ihit)**2 +  Pin(3,Ihit)**2 )
  67 C
  68       If (Abs(Pin(3,Ihit))/Pp.lt.0.0000001) go to 1000 ! Tracks paralell CPV
  69 C
  70       CorX=    2.*DxCPV/Pin(3,Ihit)
  71       Dzx = Pin(1,Ihit)*CorX            !  Incline tracks
  72       Dyx = Pin(2,Ihit)*CorX
  73 C                                       !  Charge of inline tracks
  74       CALL RNORML(R1,1)
  75       Pin(4,Ihit) =  Pin(4,Ihit)*(1.+Detr*R1)*
  76      +               SQRT(Dzx**2 + Dyx**2 + 4.*DxCPV**2)/(2.*DxCPV)
  77 C
  78 c      call hfill(113,Dzx, 0., 1.)
  79 c      call hfill(114,Dyx, 0., 1.)
  80 C
  81       Zhit1 = Xin(1,Ihit)+DzCPV
  82       Yhit1 = Xin(2,Ihit)+DyCPV         !  Initial Y coordinate
  83 C
  84       Zhit2 = Zhit1+Dzx                 !  final Y coordinate
  85       Yhit2 = Yhit1+Dyx
  86 C
  87       Iwht1 = Yhit1/CelWr               !  Wire (Y) coordinate "in"
  88       Iwht2 = Yhit2/CelWr               !  Wire (Y) coordinate "out"
  89 C
  90       If(Iwht1.eq.Iwht2)   then         !  Incline 1 wire hit
  91                            Niter      = 2
  92                            ZYEgm(1,1) =(Zhit1+Zhit2-Dzx*.57735)/2.
  93                            ZYEgm(2,1) =(Iwht1+0.5)*CelWr           !
  94                            ZYEgm(3,1) = Pin(4,Ihit)/2.
  95 C
  96                            ZYEgm(1,2) =(Zhit1+Zhit2+Dzx*.57735)/2.
  97                            ZYEgm(2,2) =(Iwht1+0.5)*CelWr           !
  98                            ZYEgm(3,2) = Pin(4,Ihit)/2.
  99                         go to 10
 100                   endif
 101 C                                       !    3 Wire tracks   !
 102 C                                     --------------------------
 103       IF(Iwht1-1.ne.Iwht2 .and.
 104      +   Iwht1+1.ne.Iwht2)     THEN
 105                                Niter = 3
 106                                Iwht3 =(Iwht1+Iwht2)/2
 107                                EGM   = Pin(4,Ihit)
 108                                Ywht1 =(Iwht1+0.5)*CelWr          !  Wire 1
 109                                Ywht2 =(Iwht2+0.5)*CelWr          !  Wire 2
 110                                Ywht3 =(Iwht3+0.5)*CelWr          !  Wire 3
 111                                Ywr13 =(Ywht1+Ywht3)/2.           !  Center 13
 112                                Ywr23 =(Ywht2+Ywht3)/2.           !  Center 23
 113 C
 114            DyW1  = Yhit1-Ywr13
 115            DyW2  = Yhit2-Ywr23
 116            Egm1  = ABS(DyW1)/(ABS(DyW1)+ABS(DyW2) + CelWr)
 117            Egm2  = ABS(DyW2)/(ABS(DyW1)+ABS(DyW2) + CelWr)
 118            Egm3  =    CelWr /(ABS(DyW1)+ABS(DyW2) + CelWr)
 119 C
 120            ZYEgm(1,1) =(Dzx*(Ywr13-Ywht1)/Dyx + Zhit1 + Zhit1)/2.
 121            ZYEgm(2,1) = Ywht1
 122            ZYEgm(3,1) = EGM*Egm1
 123 C
 124            ZYEgm(1,2) =(Dzx*(Ywr23-Ywht1)/Dyx + Zhit1 + Zhit2)/2.
 125            ZYEgm(2,2) = Ywht2
 126            ZYEgm(3,2) = EGM*Egm2
 127 C
 128            ZYEgm(1,3) = Dzx*(Ywht3-Ywht1)/Dyx + Zhit1
 129            ZYEgm(2,3) = Ywht3
 130            ZYEgm(3,3) = EGM*Egm3
 131 C
 132            go to 10
 133       ENDIF
 134 C
 135       EGM   = Pin(4,Ihit)
 136       Ywht1 =(Iwht1+0.5)*CelWr          !  Wire 1
 137       Ywht2 =(Iwht2+0.5)*CelWr          !  Wire 2
 138       Ywr12 =(Ywht1+Ywht2)/2.
 139 C
 140       DyW1  = Yhit1-Ywr12
 141       DyW2  = Yhit2-Ywr12
 142       Egm1  = ABS(DyW1)/(ABS(DyW1)+ABS(DyW2))
 143       Egm2  = ABS(DyW2)/(ABS(DyW1)+ABS(DyW2))
 144 C
 145       Niter      = 2
 146       ZYEgm(1,1) =(Zhit1+Zhit2-Dzx*Egm1)/2.
 147       ZYEgm(2,1) = Ywht1
 148       ZYEgm(3,1) = EGM*Egm1
 149 C
 150       ZYEgm(1,2) =(Zhit1+Zhit2+Dzx*Egm2)/2.
 151       ZYEgm(2,2) = Ywht2
 152       ZYEgm(3,2) = EGM*Egm2
 153 C
 154 c      CALL HFILL(110, Egm1,0.,1.)
 155 c      CALL HFILL(110, Egm2,0.,1.)
 156 C
 157 C                                       !  Finite size of ionisation region
 158   10  DO 113 Iter=1,Niter
 159       Zhit =  ZYEgm(1,Iter)             !
 160       Yhit =  ZYEgm(2,Iter)             !  Wire coordinate
 161       EGM  =  ZYEgm(3,Iter)
 162 C
 163       IF   ( Zhit.LE.0. .OR. Yhit.LE.0. ) GO TO 1000
 164 C
 165       Zcel = Zhit/CelZ
 166       Ycel = Yhit/CelY
 167       Izg  = Zcel
 168       Iyg  = Ycel
 169       IF ( Izg.GE.NcelZ.OR.Iyg.GE.NcelY ) GO TO 1000
 170 C!-
 171       Zc   = Zcel - Izg-.5
 172       Yc   = Ycel - Iyg-.5
 173 C
 174       Nz3  = (NgamZ+1)/2
 175       Ny3  = (NgamY+1)/2
 176       DO 112 JJ=1,NgamZ
 177       J=JJ-Nz3
 178       KZG=IZG+J
 179       IF (KZG.LE.0.OR.KZG.GT.NcelZ) GO TO 112
 180       ZG=FLOAT(J)-Zc
 181 C
 182       DO 111 II=1,NgamY
 183       I=II-Ny3
 184       KYG=IYG+I
 185       IF (KYG.LE.0.OR.KYG.GT.NcelY) GO TO 111
 186       YG=FLOAT(I)-Yc
 187 C
 188       CALL RNORML(R1,1)
 189       AEG = AMP(EGM,ZG,YG)
 190       AGM = Chit*AEG + R1*Cnoise
 191       IF  ( AGM.LE.0.)              GO TO 111
 192 C
 193       ICPV(KZG,KYG) = ICPV(KZG,KYG) + AGM
 194       ECPV  = ECPV  +                 AGM
 195 C
 196 C-    type *, ' KZG,KYG,AGM=',KZG,KYG,AGM
 197 C
 198       IF  (KZG.LE.NZMN) NZMN = KZG
 199       IF  (KZG.GE.NZMX) NZMX = KZG
 200       IF  (KYG.LE.NYMN) NYMN = KYG
 201       IF  (KYG.GE.NYMX) NYMX = KYG
 202  111  CONTINUE
 203  112  CONTINUE
 204  113  CONTINUE
 205       Ncpv = Ncpv + 1
 206 C
 207  1000 CONTINUE
 208 c      CALL HF1( 2,ECPV,1.)
 209 c      CALL HF1(31,ECPV,1.)
 210 C
 211       CALL COMPRS(NZMN,NcelZ,NYMN,NcelY,ICPV,NWGAMS)
 212       CALL RECONS(NWGAMS)
 213 C
 214       IF(NGAM.LE.0) RETURN
 215 C                                                !   O
 216       DO Ig=1,NGAM                               !   U
 217       Xout(1,Ig)  = XGAM(Ig) + 0.5*CelZ - DzCPV  !   T
 218       Xout(2,Ig)  = YGAM(Ig) + 0.5*CelY - DyCPV  !   P
 219       ENDDO                                      !   U
 220       Nout = NGAM                                !   T
 221 C
 222       RETURN
 223       END
 224 C=======================================================================
 225       SUBROUTINE COMPRS(NXMN,NXMX,NYMN,NYMX,IA,NWGAMS)
 226       COMMON/EVENT/ IBH(5),IPHD,IPSC,IPGS,IPGM,IPTG,I11,I12,IBUF(30035)
 227       COMMON/CONSTA/ NcelZ,NcelY,IDELT,IDELA,CelZ,CelY
 228       INTEGER*2  IPH,IA(NXMX,NYMX)
 229       DATA  IPGM,Mdim,IBUF(3) / 35, 30035, 0 /
 230 C
 231       IPGMS =IPGM
 232       DO 10 IY=NYMN,NYMX
 233       DO 10 IX=NXMN,NXMX
 234 C
 235       IPH = IA(IX,IY)
 236       IF (IPH.LT.IDELT) GO TO 10
 237 c      CALL HFILL(9,FLOAT(IPH),0.,1.)
 238 C
 239       IPGMS =IPGMS+2
 240       IF(IPGMS.ge.Mdim) GO TO 10
 241 C
 242       ICNT  =NYMX* IX   +IY
 243       IBUF(  IPGMS) = IPH
 244       IBUF(1+IPGMS) = ICNT
 245    10 CONTINUE
 246 C
 247       NWGAMS = IPGMS-IPGM
 248
 249       IBUF(1)= IPGMS
 250       RETURN
 251       END
 252 C
 253 C=======================================================================
 254       SUBROUTINE RECONS(NWGAMS)
 255       COMMON/EVENT/ IHDR(12),IBUF(30035)
 256       EQUIVALENCE(IHDR(3),NREC),(IHDR(5),NMT),(IHDR(6),IPHOD),(IHDR(7),
 257      ,IPSC),(IHDR(8),IPGS),(IHDR(9),IPGM),(IHDR(10),IPTG),(IBUF(3),NEV)
 258       COMMON/GAMMA/ NGAM,EGAM(5000),XGAM(5000),YGAM(5000),CHGAM(5000)
 259       COMMON/CLUSTER/ NCL,LENCL(5000)
 260       COMMON/WORK/ IWRK(40000)
 261       DATA NWMIN,NWMAX/2,30000/,MINPK,MINECL/10,15/,ECUT/25./
 262       DATA SCALMS/26./
 263 C
 264       NGAM=0
 265 c      CALL HFILL(1,  FLOAT(NWGAMS),  0.,1.)   !..my..!
 266 c      CALL HFILL(101,FLOAT(NWGAMS/2),0.,1.)
 267
 268       IF(NWGAMS.LT.NWMIN.OR.NWGAMS.GT.NWMAX) RETURN
 269       IPGAMS=IPGM+2
 270       ETOT=0.
 271       IPEND=IPGAMS+NWGAMS-2
 272       DO 1 I=IPGAMS,IPEND,2
 273     1 ETOT=ETOT+IBUF(I)
 274 c      CALL HFILL( 7,ETOT,0.,1.)
 275 c      CALL HFILL(32,ETOT,0.,1.)
 276       IF(ETOT.LT.ECUT)         RETURN
 277
 278       CALL CLUST(NWGAMS,IBUF(IPGAMS))
 279 c      CALL HFILL(3,FLOAT(NCL),0.,1.)  !..my..!
 280 C
 281       CALL VARDEL(NWGAMS)
 282       IF(NCL.LT.1)             RETURN
 283       ETOTCL=0.
 284       IPCL=IPGAMS
 285       DO 10 ICL=1,NCL
 286          IECL=0
 287 C
 288 c         CALL HFILL(4,FLOAT(LENCL(ICL)),0.,1.) !..my..!
 289 C
 290          DO 5 I=1,LENCL(ICL)
 291  5       IECL=IECL+IBUF(IPCL+2*I-2)
 292          IF (IECL.GE.MINECL)       ETOTCL=ETOTCL+IECL
 293 *
 294 c         IF (iecl.GE.900) THEN
 295 c            CALL hf1 (104,float(lencl(icl)),1.)
 296 c         END IF
 297 *
 298          IPCL=IPCL+LENCL(ICL)*2
 299  10   CONTINUE
 300 c      CALL HFILL( 5,ETOTCL,0.,1.)  !..my..!
 301 c      CALL HFILL(33,ETOTCL,0.,1.)  !..my..!
 302       IF(ETOTCL.LT.ECUT)       RETURN
 303 C
 304       IPCL=IPGAMS
 305       DO 30 ICL=1,NCL
 306       IECL=0
 307       DO 25 I=1,LENCL(ICL)
 308    25 IECL=IECL+IBUF(IPCL+2*I-2)
 309       IF(IECL.LT.MINECL)       GO TO 30
 310       CALL GAMS(IBUF(IPCL),LENCL(ICL),IWRK(1),MINPK)
 311    29 IF(NGAM.GE.4999)         GO TO 50
 312    30 IPCL=IPCL+LENCL(ICL)*2
 313 C
 314    50 CONTINUE
 315       CALL KILLG(SCALMS)
 316 c      CALL HFILL(6,FLOAT(NGAM),0.,1.)  !..my..!
 317       RETURN
 318       END
 319 C=======================================================================
 320       SUBROUTINE KILLG(SCALMS)
 321       COMMON /GAMMA / NGAM,EGAM(5000),XGAM(5000),YGAM(5000),CHGAM(5000)
 322       COMMON /CONSTA/ NcelZ,NcelY,IDELT,IDELA,CelZ,CelYC
 323 cc    DATA   SQDCUT / 3.50/, EFMCUT / 5./, THR0 /500./  ! N_c = 1500
 324       DATA   SQDCUT / 4.50/, EFMCUT / 5./, THR0 /600./  ! N-c = 800
 325 C
 326       NGM=0
 327       ETOT=0.
 328       DO 1 IG=1,NGAM
 329          IF(EGAM(IG).LT.THR0)                 GO TO 1
 330          NGM=NGM+1
 331          XGAM(NGM) = XGAM(IG)
 332          YGAM(NGM) = YGAM(IG)
 333          EGAM(NGM) = EGAM(IG)
 334          CHGAM(NGM)=CHGAM(IG)
 335          ETOT=ETOT + EGAM(IG)
 336 c         CALL HF1(35,EGAM(NGM),1.)
 337     1 CONTINUE
 338       NGAM=NGM
 339 c      CALL HF1(37,FLOAT(NGM),1.)
 340 C
 341       iloop = 0
 342     2 IF(NGAM.LT.2)                           RETURN
 343       SQDMIN = 1.0E+10
 344       iloop = iloop + 1
 345       DO IG1=1,NGAM-1
 346          DO IG2=IG1+1,NGAM
 347             SQDIST = (XGAM(IG1)-XGAM(IG2))**2 + (YGAM(IG1)-YGAM(IG2))**2
 348 c            IF (iloop .EQ. 1) THEN
 349 c               CALL HF1(36, SQDIST, 1.)
 350 c            END IF
 351             IF (SQDIST .LT. SQDMIN) THEN
 352                SQDMIN = SQDIST
 353                IG1M   = IG1
 354                IG2M   = IG2
 355             END IF
 356          END DO
 357       END DO
 358 C
 359       IF (SQDMIN .GT. SQDCUT) RETURN
 360 *
 361       IF (egam(ig1m) .GT. egam(ig2m)) THEN
 362          XGAM(IG1M)  =  XGAM(IG1M)
 363          YGAM(IG1M)  =  YGAM(IG1M)
 364          CHGAM(IG1M) = CHGAM(IG1M)
 365          EGAM(IG1M)  =  EGAM(IG1M)
 366       ELSE
 367          XGAM(IG1M)  =  XGAM(IG2M)
 368          YGAM(IG1M)  =  YGAM(IG2M)
 369          CHGAM(IG1M) = CHGAM(IG2M)
 370          EGAM(IG1M)  =  EGAM(IG2M)
 371       END IF
 372 c      DE          =  EGAM(IG1M)/(EGAM(IG1M)+EGAM(IG2M))
 373 c      XGAM(IG1M)  =  XGAM(IG1M)*DE +  XGAM(IG2M)*(1.-DE)
 374 c      YGAM(IG1M)  =  YGAM(IG1M)*DE +  YGAM(IG2M)*(1.-DE)
 375 c      CHGAM(IG1M) = CHGAM(IG1M)    + CHGAM(IG2M)
 376 c      EGAM(IG1M)  =  EGAM(IG1M)    +  EGAM(IG2M)
 377
 378 c      CALL hf1 (38, egam(ig1m), 1.)
 379       NGAM=NGAM-1
 380       IF(IG2M.GT.NGAM)                        GO TO 2
 381       DO IG=IG2M,NGAM
 382          XGAM(IG) = XGAM(IG+1)
 383          YGAM(IG) = YGAM(IG+1)
 384          EGAM(IG) = EGAM(IG+1)
 385          CHGAM(IG)=CHGAM(IG+1)
 386       END DO
 387       GO TO 2
 388 *
 389       END
 390 *=======================================================================
 391       SUBROUTINE ORDER(N,IAR)
 392       DIMENSION IAR(N)
 393       IF(N.LT.4) RETURN
 394       DO 4 K=4,N,2
 395       IF(IAR(K).GE.IAR(K-2)) GO TO 4
 396       IA=IAR(K)
 397       ID=IAR(K-1)
 398       IF(K.EQ.4) GO TO 2
 399       DO 1 I1=2,K-4,2
 400       I=K-I1
 401       IF(IAR(I-2).LE.IA) GO TO 3
 402       IAR(I+2)=IAR(I)
 403     1 IAR(I+1)=IAR(I-1)
 404     2 I=2
 405     3 IAR(I+2)=IAR(I)
 406       IAR(I+1)=IAR(I-1)
 407       IAR(I)=IA
 408       IAR(I-1)=ID
 409     4 CONTINUE
 410       RETURN
 411       END
 412 C=======================================================================
 413       SUBROUTINE CLUST(NW,IA)
 414       DIMENSION IA(2,NW)
 415       COMMON/CLUSTER/ NCL,LENCL(5000)
 416       COMMON/CONSTA/ NcelZ,NcelY,IDELT,IDELA,CelZ,CelY
 417       COMMON/WORK/ IWORK(40000)
 418       DATA MxCL  / 5000 /
 419       NCL=0
 420       IF(NW.LT.2) RETURN
 421       NCL=1
 422       LENCL(1)=1
 423       IF(NW.LT.4) RETURN
 424       CALL ORDER(NW,IA(1,1))
 425       NW2=NW/2
 426       NCL=0
 427       NEXT=IA(2,1)
 428       IAK=NEXT
 429       K2 =2
 430     5 DO 10 K=K2,NW2
 431       IAK1=IAK
 432       IAK =IA(2,K)
 433 C
 434       IF(IAK.EQ.IAK/NcelY*NcelY) GO TO 20
 435       IF(IAK.GT.IAK1+1)          GO TO 20
 436    10 CONTINUE
 437 C
 438    20 IBEG=NEXT
 439       NEXT=IAK
 440       IEND=IA(2,K-1)
 441       KB  = K-1 - IEND + IBEG
 442       IF(NCL.GE.MxCL)            RETURN
 443 C
 444       NCL=NCL+1
 445       LENCL(NCL)=IEND-IBEG+1
 446       IF(NCL.EQ.1.AND.K.GT.NW2)  RETURN
 447       IF(NCL.EQ.1)   THEN
 448                      K2=K+1
 449                      GO TO 5
 450                      ENDIF
 451       LAST=KB-1
 452       NCL0=NCL
 453       DO 60 ICL1=1,NCL0-1
 454       ICL=NCL0-ICL1
 455       IF(IBEG-IA(2,LAST).GT.NcelY.AND.K.GT.NW2) RETURN
 456       IF(IBEG-IA(2,LAST).GT.NcelY)  THEN
 457                                     K2=K+1
 458                                     GO TO 5
 459                                     ENDIF
 460       LENG=LENCL(ICL)
 461 C-
 462       DO 50 I=1,LENG
 463       ICUR=IA(2,LAST-I+1)
 464       IF(IEND-ICUR.LT.NcelY)        GO TO 50
 465       IF(IBEG-ICUR.GT.NcelY)        GO TO 60
 466       IF(ICL.EQ.NCL-1)           GO TO 40
 467       IF(LENG.GT.768)            GO TO 60
 468       CALL UCOPY(IA(1,LAST+1-LENG),IWORK(1),LENG*2)
 469       CALL UCOPY(IA(1,LAST+1),IA(1,LAST+1-LENG),(KB-1-LAST)*2)
 470       CALL UCOPY(IWORK(1),IA(1,KB-LENG),LENG*2)
 471       DO 30 J=ICL,NCL-2
 472    30 LENCL(J)=LENCL(J+1)
 473    40 NCL=NCL-1
 474       KB = KB - LENG
 475       LENCL(NCL)=K-KB
 476       GO TO 60
 477    50 CONTINUE
 478    60 LAST=LAST-LENG
 479       IF(K.LE.NW2)   THEN
 480                      K2=K+1
 481                      GO TO 5
 482                      ENDIF
 483       RETURN
 484       END
 485 C=======================================================================
 486       SUBROUTINE GAMS(IGAMS,NADC,IWRK,MINPK)
 487 C  15.MAR.1983.
 488       COMMON/GAMMA/ NGAM,EGAM(5000),XGAM(5000),YGAM(5000),CHGAM(5000)
 489       COMMON/CONSTA/ NcelZ,NcelY,IDELT,IDELA,CelZ,CelY
 490       DIMENSION IPNPK(50),EPK(50),XPK(50),YPK(50),IGMPK(2,50)
 491       DIMENSION IGAMS(2,NADC),IWRK(NADC,50)
 492       DATA MxPK                       / 50/,CHISQ1/30./,CHISQ2/3./
 493 C
 494       IDELTA = IDELA
 495 C
 496       CALL ORDER(NADC*2,IGAMS(1,1))
 497       NGAM0=NGAM
 498 C
 499 C   PEAKS SEARCH
 500 C
 501       NPK=0
 502       DO 50 IC=1,NADC
 503       IAC=IGAMS(1,IC)
 504       IF(IAC.LT.MINPK)       GO TO 50
 505       IXY=IGAMS(2,IC)
 506       IXYC=IXY
 507       IYC=IXYC-IXYC/NcelY*NcelY
 508       IF(NADC.LT.2)          GO TO 40
 509       IF(IC.EQ.NADC)         GO TO 20
 510       DO 10 IN=IC+1,NADC
 511       IXY=IGAMS(2,IN)
 512       IF(IXY-IXYC.GT.NcelY+1)    GO TO 20
 513       IF(IABS(IXY-IXY/NcelY*NcelY-IYC).GT.1)  GO TO 10
 514       IF(IGAMS(1,IN).GE.IAC) GO TO 50
 515    10 CONTINUE
 516    20 IF(IC.LT.2)            GO TO 40
 517       DO 30 IN1=1,IC-1
 518       IN=IC-IN1
 519       IXY=IGAMS(2,IN)
 520       IF(IXYC-IXY.GT.NcelY+1)     GO TO 40
 521       IF(IABS(IXY-IXY/NcelY*NcelY-IYC).GT.1)  GO TO 30
 522       IF(IGAMS(1,IN).GT.IAC) GO TO 50
 523    30 CONTINUE
 524    40 NPK=NPK+1
 525       IPNPK(NPK)=IC
 526 C!!-  IF(NPK.EQ.  10.OR.NPK.GE.NcelY*NcelY/NADC-3) GO TO 60
 527       IF(NPK.EQ.MxPK.OR.NPK.GE.NcelZ*NcelY/NADC-3) GO TO 60
 528    50 CONTINUE
 529       IF(NPK.EQ.0)           RETURN
 530 C
 531    60 CONTINUE
 532 C
 533 C   GAMMA SEARCH FOR ONE PEAK
 534 C
 535   101 IF(NPK.GT.1)           GO TO 102
 536       IF(NGAM.GE.4999)       RETURN
 537       NGAM=NGAM+1
 538       CHISQ=CHISQ2
 539       CALL GAM (NADC,IGAMS(1,1),CHISQ,EGAM(NGAM),XGAM(NGAM),YGAM(NGAM),
 540      ,                          EGAM(NGAM+1),XGAM(NGAM+1),YGAM(NGAM+1))
 541       CHGAM(NGAM)=CHISQ
 542       IF(EGAM(NGAM+1).EQ.0..OR.NGAM.EQ.4999)  GO TO 500
 543       NGAM=NGAM+1
 544       CHGAM(NGAM)=CHISQ
 545       GO TO 500
 546 C
 547 C   GAMMA SEARCH FOR SEVERAL PEAKS
 548 C   FIRST STEP AND STEP 1-BIS.
 549 C
 550   102 IF(NGAM.GE.4999)         RETURN
 551       DO 170 IBIS=1,2
 552       DO 105 I=1,NADC
 553   105 IWRK(I,NPK+3)=0
 554 C
 555       DO 160 IPK=1,NPK
 556       IC=IPNPK(IPK)
 557       E=IGAMS(1,IC)
 558       IF(IBIS.NE.1) E=E*FLOAT(IWRK(IC,IPK))/IWRK(IC,NPK+1)
 559       IXYPK=IGAMS(2,IC)
 560       IXPK=IXYPK/NcelY
 561       IYPK=IXYPK-IXPK*NcelY
 562       EPK(IPK)=E
 563       XPK(IPK)=E*IXPK
 564       YPK(IPK)=E*IYPK
 565       IF(IC.GE.NADC)         GO TO 120
 566       DO 110 IN=IC+1,NADC
 567       IXY=IGAMS(2,IN)
 568       IF(IXY-IXYPK.GT.NcelY+1)    GO TO 120
 569       IX=IXY/NcelY
 570       IY=IXY-IX*NcelY
 571       IF(IABS(IY-IYPK).GT.1) GO TO 110
 572       E=IGAMS(1,IN)
 573       IF(IBIS.NE.1) E=E*FLOAT(IWRK(IN,IPK))/IWRK(IN,NPK+1)
 574       EPK(IPK)=EPK(IPK)+E
 575       XPK(IPK)=XPK(IPK)+E*IX
 576       YPK(IPK)=YPK(IPK)+E*IY
 577   110 CONTINUE
 578   120 IF(IC.LT.2)            GO TO 140
 579       DO 130 IN1=1,IC-1
 580       IN=IC-IN1
 581       IXY=IGAMS(2,IN)
 582       IF(IXYPK-IXY.GT.NcelY+1)    GO TO 140
 583       IX=IXY/NcelY
 584       IY=IXY-IX*NcelY
 585       IF(IABS(IY-IYPK).GT.1) GO TO 130
 586       E=IGAMS(1,IN)
 587       IF(IBIS.NE.1) E=E*FLOAT(IWRK(IN,IPK))/IWRK(IN,NPK+1)
 588       EPK(IPK)=EPK(IPK)+E
 589       XPK(IPK)=XPK(IPK)+E*IX
 590       YPK(IPK)=YPK(IPK)+E*IY
 591   130 CONTINUE
 592   140 CONTINUE
 593       XPK(IPK)=XY(XPK(IPK)/EPK(IPK))
 594       YPK(IPK)=XY(YPK(IPK)/EPK(IPK))
 595       DO 150 I=1,NADC
 596       IXY=IGAMS(2,I)
 597       IX=IXY/NcelY
 598       IY=IXY-IX*NcelY
 599       DX=IX-XPK(IPK)
 600       DY=IY-YPK(IPK)
 601       IWK=0
 602       IF(DX*DX+DY*DY.LT.8.) IWK = AMP(EPK(IPK),DX,DY)+.5
 603       IWRK(I,IPK)=IWK
 604       IWRK(I,NPK+3)=IWRK(I,NPK+3)+IWK
 605   150 CONTINUE
 606   160 CONTINUE
 607       DO 165 I=1,NADC
 608       IWK=IWRK(I,NPK+3)
 609       IF(IWK.LE.0) IWK=1
 610   165 IWRK(I,NPK+1)=IWK
 611   170 CONTINUE
 612 C
 613       DO 200 IPK=1,NPK
 614       LENG=0
 615       DO 190 I=1,NADC
 616       IF(IWRK(I,NPK+3).LE.0) GO TO 190
 617       IE=IGAMS(1,I)*FLOAT(IWRK(I,IPK))/IWRK(I,NPK+3)+.5
 618       IF(IE.LE.IDELTA)       GO TO 190
 619       LENG=LENG+1
 620       IWRK(2*LENG-1,NPK+1)=IE
 621       IWRK(2*LENG,NPK+1)=IGAMS(2,I)
 622   190 CONTINUE
 623       IF(NGAM.GE.4999)       GO TO 500
 624       IF(LENG.EQ.0)          GO TO 200
 625       NGAM=NGAM+1
 626       CHISQ=CHISQ1
 627       CALL GAM (LENG,IWRK(1,NPK+1),CHISQ,EGAM(NGAM),XGAM(NGAM),
 628      ,       YGAM(NGAM),EGAM(NGAM+1),XGAM(NGAM+1),YGAM(NGAM+1))
 629       CHGAM(NGAM)=CHISQ
 630       IGMPK(1,IPK)=NGAM
 631       IGMPK(2,IPK)=NGAM
 632       IF(EGAM(NGAM+1).EQ.0..OR.NGAM.EQ.4999)   GO TO 200
 633       NGAM=NGAM+1
 634       CHGAM(NGAM)=CHISQ
 635       IGMPK(2,IPK)=NGAM
 636   200 CONTINUE
 637 C
 638 C   SECOND STEP.
 639 C
 640       DO 210 I=1,NADC
 641   210 IWRK(I,NPK+3)=0
 642       DO 230 IPK=1,NPK
 643       DO 230 I=1,NADC
 644       IWRK(I,IPK)=0
 645       DO 220 IG=IGMPK(1,IPK),IGMPK(2,IPK)
 646       IXY=IGAMS(2,I)
 647       DX=IXY/NcelY-XGAM(IG)
 648       DY=IXY-IXY/NcelY*NcelY-YGAM(IG)
 649       IF(DX*DX+DY*DY.GT.8.)  GO TO 220
 650       IWK = AMP(EGAM(IG),DX,DY)+.5
 651       IWRK(I,IPK)  =IWRK(I,IPK)  +IWK
 652       IWRK(I,NPK+3)=IWRK(I,NPK+3)+IWK
 653   220 CONTINUE
 654   230 CONTINUE
 655 C
 656       NGAM=NGAM0
 657       DO 300 IPK=1,NPK
 658       LENG=0
 659       DO 290 I=1,NADC
 660       IF(IWRK(I,NPK+3).LE.0) GO TO 290
 661       IE=IGAMS(1,I)*FLOAT(IWRK(I,IPK))/IWRK(I,NPK+3)+.5
 662       IF(IE.LE.IDELTA)       GO TO 290
 663       LENG=LENG+1
 664       IWRK(2*LENG-1,NPK+1)=IE
 665       IWRK(2*LENG,NPK+1)=IGAMS(2,I)
 666   290 CONTINUE
 667       IF(NGAM.GE.4999)       GO TO 500
 668       IF(LENG.EQ.0)          GO TO 300
 669       NGAM=NGAM+1
 670       CHISQ=CHISQ2
 671       CALL GAM (LENG,IWRK(1,NPK+1),CHISQ,EGAM(NGAM),XGAM(NGAM),
 672      ,       YGAM(NGAM),EGAM(NGAM+1),XGAM(NGAM+1),YGAM(NGAM+1))
 673       CHGAM(NGAM)=CHISQ
 674       IF(EGAM(NGAM+1).EQ.0..OR.NGAM.EQ.4999)   GO TO 300
 675       NGAM=NGAM+1
 676       CHGAM(NGAM)=CHISQ
 677   300 CONTINUE
 678   500 DO 600 IG=NGAM0+1,NGAM
 679       XGAM(IG)= XGAM(IG)*CelZ
 680   600 YGAM(IG)= YGAM(IG)*CelY
 681       RETURN
 682       END
 683 C=======================================================================
 684       SUBROUTINE GAM(NADC,IADC,CHISQ,E1,X1,Y1,E2,X2,Y2)
 685 C   28 OCT.1981.     ONE GAMMA FIT PROGRAM.
 686 CORRECTED 28 MAR.1983.
 687       DIMENSION IADC(2,NADC)
 688       DATA STPMIN/.005/,CONST/1.5/,ST0/.00005/,CHMIN/1./
 689       COMMON/CONSTA/ NcelZ,NcelY,IDELT,IDELA,CelZ,CelY
 690 C
 691 C     WRITE(5,1)
 692 C   1 FORMAT(' *** GAM ***')
 693       IF(NADC.LE.0)                   RETURN
 694       DOF=NADC-2
 695       IF(DOF.LT.1.)                   DOF=1.
 696       CHSTOP=CHMIN*DOF
 697       CSQCUT=CHISQ
 698       E1=0.
 699       XX=0.
 700       YY=0.
 701       E2=0.
 702       X2=0.
 703       Y2=0.
 704       DO 5 I=1,NADC
 705       A  = IADC(1,I)
 706       E1 = E1 + A
 707       IXY=IADC(2,I)
 708       IX=IXY/NcelY
 709       IY=IXY-IX*NcelY
 710       XX = XX + A*IX
 711       YY = YY + A*IY
 712     5 CONTINUE
 713       XC = XY(XX/E1)
 714       YC = XY(YY/E1)
 715       ST = ST0
 716       CHISQ=1.E20
 717       GR = 1.
 718       GRX=0.
 719       GRY=0.
 720       DO 100 ITER=1,50
 721       CHISQC=0.
 722       GRXC  =0.
 723       GRYC  =0.
 724       DO 10 I=1,NADC
 725       IXY=IADC(2,I)
 726       DX = XC - IXY/NcelY
 727       DY = YC - (IXY-IXY/NcelY*NcelY)
 728       CALL AG(E1,DX,DY,A,GX,GY)
 729       D  = CONST * A * (1.-A/E1)
 730       IF(D.LT.0.) D = 0.
 731       D  =     9.
 732       EX = IADC(1,I)
 733       DA = A - EX
 734       CHISQC=CHISQC+DA**2/D
 735       DD = DA/D
 736       DD = DD*(2.-DD*CONST*(1.-2.*A/E1))
 737       GRXC = GRXC + GX*DD
 738       GRYC = GRYC + GY*DD
 739    10 CONTINUE
 740       GRC  = SQRT(GRXC**2+GRYC**2)
 741       IF(GRC.LT.1.E-10) GRC=1.E-10
 742       SC = 1.+CHISQC/CHISQ
 743       ST = ST/SC
 744       IF(CHISQC.GT.CHISQ)              GO TO 20
 745       COSI=(GRX*GRXC+GRY*GRYC)/GR/GRC
 746       ST = ST*SC/(1.4-COSI)
 747       X1 = XC
 748       Y1 = YC
 749       GRX=GRXC
 750       GRY=GRYC
 751       GR = GRC
 752       CHISQ=CHISQC
 753       IF(CHISQ.LT.CHSTOP)              GO TO 101
 754    20 STEP=ST*GR
 755       IF(STEP.LT.STPMIN)               GO TO 101
 756       IF(STEP.GT.1.)                   ST = 1./GR
 757       XC = X1 - ST*GRX
 758       YC = Y1 - ST*GRY
 759   100 CONTINUE
 760   101 CHISQ=CHISQ/DOF
 761       IF(CHISQ.GT.CSQCUT) CALL TWOGAM(NADC,IADC,CHISQ,E1,X1,Y1,E2,X2,Y2)
 762       RETURN
 763       END
 764 C=======================================================================
 765       SUBROUTINE TWOGAM(NADC,IADC,CHISQ,E1,X1,Y1,E2,X2,Y2)
 766 C   28 OCT.1981.     TWO GAMMA FIT PROGRAM.
 767 CORRECTED 28 MAR.83.
 768       COMMON/CONSTA/ NcelZ,NcelY,IDELT,IDELA,CelZ,CelY
 769       DIMENSION IADC(2,NADC)
 770       DATA STPMIN/.005/,CONST/1.5/,DELCH/6./,EMIN/20./,ST0/.00005/
 771       DATA CHMIN/1./
 772 C
 773 C     WRITE(5,7)
 774 C   7 FORMAT(1X,'*** TWOGAM ***')
 775       IF(NADC.LE.3)                   RETURN
 776       DOF=NADC-5
 777       IF(DOF.LT.1.)                   DOF=1.
 778       CHSTOP=CHMIN*DOF
 779 C   START POINT CHOOSING.
 780       E=0.
 781       XX=0.
 782       YY=0.
 783       DO 1 I=1,NADC
 784       A=IADC(1,I)
 785       IXY=IADC(2,I)
 786       IX=IXY/NcelY
 787       IY=IXY-IX*NcelY
 788       XX=XX + A*IX
 789       YY=YY + A*IY
 790       E = E + A
 791     1 CONTINUE
 792       XX=XY(XX/E)
 793       YY=XY(YY/E)
 794       EXX=0.
 795       EYY=0.
 796       EXY=0.
 797       DO 2 I=1,NADC
 798       A=IADC(1,I)
 799       IXY=IADC(2,I)
 800       IX=IXY/NcelY
 801       IY=IXY-IX*NcelY
 802       DX=IX-XX
 803       DY=IY-YY
 804       EXX=EXX + A*DX**2
 805       EYY=EYY + A*DY**2
 806       EXY=EXY + A*DX*DY
 807     2 CONTINUE
 808       SQR=SQRT(4.*EXY**2+(EXX-EYY)**2)
 809       EUU=(EXX+EYY+SQR)/2.
 810       COS2FI=1.
 811       IF(SQR.GT.1.E-10) COS2FI=(EXX-EYY)/SQR
 812       COSFI=SQRT((1.+COS2FI)/2.)
 813       SINFI=SQRT((1.-COS2FI)/2.)
 814       IF(EXY.LT.0.) SINFI=-SINFI
 815       EU3=0.
 816       DO 3 I=1,NADC
 817       A=IADC(1,I)
 818       IXY=IADC(2,I)
 819       IX=IXY/NcelY
 820       IY=IXY-IX*NcelY
 821       DX=IX-XX
 822       DY=IY-YY
 823       DU=DX*COSFI+DY*SINFI
 824     3 EU3=EU3+A*DU**3
 825       C=E*EU3**2/EUU**3/2.
 826       SIGN=1.
 827       IF(EU3.LT.0.) SIGN=-1.
 828       R=1.+C+SIGN*SQRT(2.*C+C**2)
 829       IF(ABS(R-1.).GT..1) U=EU3/EUU*(R+1.)/(R-1.)
 830       IF(ABS(R-1.).LE..1) U=SQRT(SQR/E/R)*(1.+R)
 831       E2C=E/(1.+R)
 832       E1C=E-E2C
 833       U1=-U/(1.+R)
 834       U2=U+U1
 835       X1C=XX+U1*COSFI
 836       Y1C=YY+U1*SINFI
 837       X2C=XX+U2*COSFI
 838       Y2C=YY+U2*SINFI
 839 C
 840 C   START OF ITERATIONS.
 841 C
 842       ST = ST0
 843       CH = 1.E20
 844       GRX1 = 0.
 845       GRY1 = 0.
 846       GRX2 = 0.
 847       GRY2 = 0.
 848       GRE  = 0.
 849       GR   = 1.
 850       DO 100 ITER=1,50
 851       CHISQC= 0.
 852       GRX1C = 0.
 853       GRY1C = 0.
 854       GRX2C = 0.
 855       GRY2C = 0.
 856       GREC  = 0.
 857       DO 10 I=1,NADC
 858       IXY=IADC(2,I)
 859       IX=IXY/NcelY
 860       IY=IXY-IX*NcelY
 861       DX1 = X1C-IX
 862       DY1 = Y1C-IY
 863       CALL AG(E1C,DX1,DY1,A1,GX1,GY1)
 864       DX2 = X2C-IX
 865       DY2 = Y2C-IY
 866       CALL AG(E2C,DX2,DY2,A2,GX2,GY2)
 867       EX = IADC(1,I)
 868       A  = A1 + A2
 869       D  = CONST * A * (1.-A/E)
 870       IF(D.LT.0.) D = 0.
 871       D  =     9.
 872       DA = A - EX
 873       CHISQC=CHISQC+DA**2/D
 874       DD = DA/D
 875       DD = DD*(2.-DD*CONST*(1.-2.*A/E))
 876       GRX1C = GRX1C + GX1*DD
 877       GRY1C = GRY1C + GY1*DD
 878       GRX2C = GRX2C + GX2*DD
 879       GRY2C = GRY2C + GY2*DD
 880       GREC  = GREC  + (A1/E1C-A2/E2C)*E*DD
 881    10 CONTINUE
 882       GRC=SQRT(GRX1C**2+GRY1C**2+GRX2C**2+GRY2C**2+GREC**2)
 883       IF(GRC.LT.1.E-10) GRC=1.E-10
 884       SC = 1.+CHISQC/CH
 885       ST = ST/SC
 886       IF(CHISQC.GT.CH)              GO TO 20
 887       COSI=(GRX1*GRX1C+GRY1*GRY1C+GRX2*GRX2C+GRY2*GRY2C+GRE*GREC)/GR/GRC
 888       ST = ST*SC/(1.4-COSI)
 889       EE1 = E1C
 890       XX1 = X1C
 891       YY1 = Y1C
 892       EE2 = E2C
 893       XX2 = X2C
 894       YY2 = Y2C
 895       CH  = CHISQC
 896       IF(CH.LT.CHSTOP)               GO TO 101
 897       GRX1= GRX1C
 898       GRY1= GRY1C
 899       GRX2= GRX2C
 900       GRY2= GRY2C
 901       GRE = GREC
 902       GR  = GRC
 903    20 STEP= ST*GR
 904       IF(STEP.LT.STPMIN)            GO TO 101
 905       DE  = ST  * GRE * E
 906       E1C = EE1 - DE
 907       E2C = EE2 + DE
 908       IF(E1C.GT.EMIN.AND.E2C.GT.EMIN) GO TO 25
 909       ST  = ST / 2.
 910       GO TO 20
 911    25 X1C = XX1 - ST*GRX1
 912       Y1C = YY1 - ST*GRY1
 913       X2C = XX2 - ST*GRX2
 914       Y2C = YY2 - ST*GRY2
 915   100 CONTINUE
 916   101 IF(CH.GE.CHISQ*(NADC-2)-DELCH)         RETURN
 917       CHISQ=CH/DOF
 918       E1 = EE1
 919       X1 = XX1
 920       Y1 = YY1
 921       E2 = EE2
 922       X2 = XX2
 923       Y2 = YY2
 924       RETURN
 925       END
 926 C=======================================================================
 927       FUNCTION XY(XI)
 928 C   MADE FROM 25 GEV ELECTRON CALIBRATION WITHOUT HODOSCOPE.  RUN APR.1981.
 929       I=XI
 930       X=XI-I-.5
 931       XI2=X*X
 932       XY=X*(.22466+XI2*(3.74014+XI2*(-25.88467+
 933      + XI2*(166.90556-XI2*294.35077))))+I+.5
 934       RETURN
 935       END
 936 C
 937 C=======================================================================
 938       SUBROUTINE  AG(Ei,Xi,Yi,Ai,GXi,GYi)
 939       REAL*4         Ei,Xi,Yi,Ai,GXi,GYi  !...Interface...!
 940       REAL*8         E ,X ,Y ,A ,GX ,GY   !..Calculation..!
 941 *  ..........................................................
 942 *  ...It uses Lednev's Amplitude calculation function.Fcml..
 943 *  ...To Calculate an Ammplitude (A) and Garadient (GX,GY)...
 944 *  ...i use REAL*4 for inteface purpose only ................
 945 *  ...calculations in REAL*8.................................
 946 *  ...................typed by Sadovsky.A.S..................
 947       REAL*8 Fcml,Dx,Dy,D
 948       INTEGER i
 949       COMMON/CONSTA/ NcelZ,NcelY,IDELT,IDELA,CelZ,CelY
 950 C
 951       Dx= CelZ/2.
 952       Dy= CelY/2.
 953 C
 954       X=Xi*CelZ  !....REAL*4....!
 955       Y=Yi*CelY  !.....to.......!
 956       E=Ei       !..............!
 957       A=Ai       !....REAL*8....!
 958       GX=GXi     !..Convertion..!
 959       GY=GYi     !..............!
 960
 961       A = Fcml(x+dx,y+dy) - Fcml(x+dx,y-dy) -
 962      +    Fcml(x-dx,y+dy) + Fcml(x-dx,y-dy)
 963       Ai= A*E
 964 C
 965       GX= GradX(x+dx,y+dy) - GradX(x+dx,y-dy) -
 966      +    GradX(x-dx,y+dy) + GradX(x-dx,y-dy)
 967       GXi=GX*E*E
 968 C
 969       GY= GradY(x+dx,y+dy) - GradY(x+dx,y-dy) -
 970      +    GradY(x-dx,y+dy) + GradY(x-dx,y-dy)
 971       GYi=GY*E*E
 972 C
 973       RETURN
 974       END
 975 C=======================================================================
 976       FUNCTION Fcml(X,Y)
 977 C                     *..Cumuliative function calculation in REAL*8..*
 978 C
 979       REAL*8   Fcml,X,Y
 980       REAL*8    A,b,Fff
 981       DATA      A,b / 1.0, 1.082 /
 982 C-                                !...Cumuliative function calculation.......!
 983 C
 984 C                                 !  Kumulyativnaya funkciya ne pravil'naya, !
 985 C                                    no poluchaemoe pri etom enrgovydelenie
 986 C                                          v yachejke PRAVIL'NOE !!!
 987 C
 988       Fff  = A*DATAN(x*y/(b*DSQRT(b*b+x*x+y*y)))
 989       Fcml = Fff/6.2831853071796D0
 990       RETURN
 991       END
 992 C=======================================================================
 993       FUNCTION GradX(X,Y)
 994 C                      *..Cumuliative function Gradient_X in REAL*8..*
 995 C
 996       REAL*8   GradX,X,Y
 997       REAL*8   a,b,skv,s,Gradient
 998       DATA     a,b / 1.0, 1.082 /
 999 C
1000       skv      = b*b + x*x + y*y
1001       Gradient = a*y*(1.-x*x/skv)*b*DSQRT(skv)/(b*b*skv +x*x*y*y)
1002       GradX    = Gradient/6.2831853071796d0
1003       RETURN
1004       END
1005 C=======================================================================
1006       FUNCTION GradY(X,Y)
1007 C                      *..Cumuliative function Gradient_Y in REAL*8..*
1008 C
1009 C-    COMMON/CONSTS/ BEAM,ANOR,ANORT,DELTA,CELL,DSST,NCNX,NCNY,STCO
1010       REAL*8   GradY,X,Y
1011       REAL*8   a,b,skv,s,Gradient
1012       DATA     a,b / 1.0, 1.082 /
1013 C
1014       skv      = b*b + x*x + y*y
1015       Gradient = a*x*(1.-y*y/skv)*b*DSQRT(skv)/(b*b*skv +x*x*y*y)
1016       GradY    = Gradient/6.2831853071796d0
1017       RETURN
1018       END
1019 C=======================================================================
1020       FUNCTION AMP(Ei,Xi,Yi)
1021 C!-
1022 C     *..Amplitude calculation in REAL*4 using cumuliative function..*
1023 C     *........to calculate amplitude only in gams cell..............*
1024 C
1025       COMMON/CONSTA/ NcelZ,NcelY,IDELT,IDELA,CelZ,CelY
1026       REAL*4   AMP,Ei,Xi,Yi,Fcm
1027       REAL*8   Fcml,E,X ,Y, Dx,Dy,D
1028       INTEGER i
1029 C
1030       Dx= CelZ/2.
1031       Dy= CelY/2.
1032       x = Xi*CelZ
1033       y = Yi*CelY
1034 C
1035       E = Ei
1036       A = Fcml(x+dx,y+dy) - Fcml(x+dx,y-dy) -
1037      +    Fcml(x-dx,y+dy) + Fcml(x-dx,y-dy)
1038       AMP = A*E
1039       RETURN
1040       END
1041
1042 C=======================================================================
1043       SUBROUTINE VARDEL(NWGAMS)
1044       COMMON/DELTA/ IDLT(40000)
1045       COMMON/CLUSTER/ NCL,LENCL(5000)
1046       COMMON/EVENT/ IHDR(12),IBUF(30035)
1047       EQUIVALENCE(IHDR(3),NREC),(IHDR(5),NMT),(IHDR(6),IPHOD),(IHDR(7),
1048      ,IPSC),(IHDR(8),IPGS),(IHDR(9),IPGM),(IHDR(10),IPTG),(IBUF(3),NEV)
1049       DATA MINDLT/2/,NTRIG/1000/,LENMAX/6/
1050 C
1051       IF(NWGAMS.LE.2.OR.NCL.LT.1)      RETURN
1052       ITRIG=ITRIG+1
1053       IF(ITRIG.LE.NTRIG)               GO TO 20
1054       ITRIG=0
1055       DO 10 I=1,1
1056       IF(IDLT(I).GT.MINDLT)            IDLT(I)=IDLT(I)-1
1057    10 CONTINUE
1058    20 IPCL=IPGM+2
1059       DO 50 ICL=1,NCL
1060       IF(LENCL(ICL).NE.1)              GO TO 30
1061       ICNT=IBUF(IPCL+1)+1
1062       IDLT(ICNT)=IDLT(ICNT)+1
1063       GO TO 50
1064    30 IF(LENCL(ICL).GT.LENMAX)         GO TO 50
1065       IPEND=IPCL+LENCL(ICL)*2-2
1066       DO 40 I=IPCL,IPEND,2
1067       ICNT=IBUF(I+1)+1
1068       IF(IDLT(ICNT).GT.MINDLT)         IDLT(ICNT)=IDLT(ICNT)-1
1069    40 CONTINUE
1070    50 IPCL=IPCL+LENCL(ICL)*2
1071       RETURN
1072       END
1073 C=======================================================================