C   28/03/97 703282257  MEMBER NAME  JEOSUM2  (JADE56.S)    FORTRAN
C   26/06/78 704172134  MEMBER NAME  @CERFITH (S)           FORTRAN
C***** CERENKOV LIGHT DETECTION EFFICIENCY FIT *****
C     S.YAMADA   15-12-77
C     RE-CODED ON 15-12-77 TO RECOVER LOST VERSION
C     UPDATED  2-10-85 TO SF6, J.OLSSON
C     FORMATS MODIFIED, 5.1.1986    J.OLSSON
C  2 COMPONENT VERSION 30.1.1986    J.OLSSON
C  A4-A7 IN CC READ IN, VARIABLE
C
C         GEP VERSION
C
C     READ CONVERTED PHOTOELECTRON NUMBERS FROM THE TABLE
C     FOR DIFFERENT DIRECTION AND POSITIONS OF THE ORIGINS.
C---- FIT THE NUMBERS FOR FIXED DIRECTION FIRST AND THEN THE FIT
C     PARAMETERS ARE FITTED AS A FUNCTION OF THE DIRECTION.
C
C   SEE DESCRIPTION IN JADE NOTE 20, BY S.YAMADA
C
C
      REAL*8 FLAMDA
      DIMENSION AVN(20),ERR(20), TRY(20),RESLT(50,13)
      DIMENSION X(20),YFIT(21), A(8),DELTAA(8),SIGMAA(8)
      DATA AVN/20*0./,ERR/20*0./,TRY/19*100.,60./,DELTAA/8*1.0/
      DATA X/0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5,
     1  13.5,14.5,15.5,16.5,17.5,18.5,19.3/
      DATA MODE/1/, ND/0/
      DIMENSION SD(20,50), SER(20,50), BC(8,4)
      LOGICAL*1 TEXT
      DIMENSION TEXT(72)
C
      COMMON /COSTB/ COST,MODEFT,AESS,BASI,CC,CC1,AEXPA,A4,A5,A6,A7
     $ ,PART1,PART2,PART3
C
      DATA IENDCP/0/
C
C---- IENDCP=1  FOR THE ENDCAP, IENDCP=0  FOR THE BARREL PART.
C
      MODEFT = 1
      NTERMS = 4
      IPRO = 0
      NPTS = 20
      IF(IENDCP.EQ.1) NPTS = 12
      IF(IENDCP.EQ.1) X(12) = 11.1
C
C PLOT INITIALIZATION
C
      CALL ACORE(10000)
      DO 3705  I = 1,21
      CALL HBOOK1( I,' NR OF PHOTONS VS. RAD. LENGTH$',20,0.,20.)
      CALL HBOOK1(30+ I,' NR OF PHOTONS FITTED, SINGLE$',20,0.,20.)
      CALL HBOOK1(60+ I,' NR OF PHOTONS FITTED, COMBINED$',20,0.,20.)
3705  CONTINUE
      DO 3706  I = 1,4
      CALL HBOOK1(100+ I,' FIT PARAMETER VS. COSTHETA$',21,-1.05,1.05)
      CALL HBOOK1(110+ I,' FIT PARAMETER, GLOBAL FIT$',21,-1.05,1.05)
3706  CONTINUE
      CALL HBOOK1(121,'CHISQ VS. COSTHETA$',21,-1.05,1.05)
      CALL HBOOK1(151,' CC1(COST)VS COSTHETA$',21,-1.05,1.05)
      CALL HBOOK1(152,' CC=10**CC1 VS COSTHETA$',21,-1.05,1.05)
      CALL HBOOK1(153,' PART1 VS COSTHETA$',21,-1.05,1.05)
      CALL HBOOK1(154,' PART2 VS COSTHETA$',21,-1.05,1.05)
      CALL HBOOK1(155,' PART3 VS COSTHETA$',21,-1.05,1.05)
      CALL HBOOK1(161,' CC1(COST)VS COSTHETA,GLOBAL$',21,-1.05,1.05)
      CALL HBOOK1(162,' CC=10**CC1 VS COSTHETA,GLOBAL$',21,-1.05,1.05)
      CALL HBOOK1(163,' PART1 VS COSTHETA,GLOBAL$',21,-1.05,1.05)
      CALL HBOOK1(164,' PART2 VS COSTHETA,GLOBAL$',21,-1.05,1.05)
      CALL HBOOK1(165,' PART3 VS COSTHETA,GLOBAL$',21,-1.05,1.05)
C
C
7608  FORMAT(72A1)
7508  FORMAT(' ',72A1)
      READ 7608, TEXT
      WRITE(6,7508) TEXT
C
C---- READ THE MONTE CARLO RESULTS.
C---- AESS,BASI,AEXPA,A4-7 ARE PARAMETERS USED IN FUNCTN
C
      READ(5,501) AESS,BASI,AEXPA,A4,A5,A6,A7
  501 FORMAT(3F10.2,4F10.4)
      WRITE(6,6610) AESS,BASI,AEXPA,A4,A5,A6,A7
6610  FORMAT('  FUNCTN PAR.: AESS BASI AEXPA ',3F10.2,' A4-7 ',4F10.4)
C
C  LOOP OVER VARIOUS COST STARTING HERE   ###########################
C
 1    CONTINUE
      READ(5,500,END=100) COST
      READ(5,502,END=100) (AVN(K),K=1,10)
      READ(5,502,END=100) (AVN(K),K=11,NPTS)
  500 FORMAT(F6.3)
  502 FORMAT(10F7.2)
C  COMPENSATE FOR INCOMPLETE STEP AT END
C     AVN(NPTS) = AVN(NPTS)*1.7
C
C---- CALCULATE THE STATISTICAL ERRORS
C---- FORMULA ASSUMES AVN IS MEAN OF TRY ATTEMPTS, ERR IS THEN
C---- ERROR ON MEAN:  SIGMA**2/TRY, WITH SIGMA = ROOT(AVN)
C
        DO 2 K=1,NPTS
    2   ERR(K) = SQRT(AVN(K)/TRY(K))*0.5
C----   THE FACTOR 0.5 IS DUE TO THE CONVERSION EFFICIENCY MAX.(0.25).
C
C  SINCE AVN IS NOT A MEAN, BUT THE SUM OF TRY ATTEMPTS, THE CORRECTER
C  ERROR IS SIMPLY SQRT(SD)
C
C       DO 2 K=1,NPTS
C   2   ERR(K) = SQRT(AVN(K))
C
C---- STORE THE SOURCE DATA FOR FUTURE USE.
      ND = ND+1
      IF(ND.GT.50) GO TO 33
        DO 3 K=1,NPTS
        SD(K,ND) = AVN(K)
    3   SER(K,ND) = ERR(K)
33     CONTINUE
C
C
C PLOT INPUT
C
C     CALL HISTOG(20,SD(1,ND),X,0.,1.,0,4,DUMMY)
C     SUBROUTINE HISTOG (NBIN,CINCON,THEORY,BINLOW,BINSIZ,KND,NL,HISTIT)
C    NBIN = NO. OF BINS
C    CINCON = EXPERIMENTAL DISTRIBUTION
C    THEORY=THEORETICAL DISTRIBUTION (SUPERIMPOSED IF REQUIRED)
C    BINLOW = LOW EDGE OF FIRST BIN
C    BINSIZ = SIZE OF BIN
C    KND=0 IF YOU DON'T WISH TO SUPPLY A THEORETICAL CURVE
C    KND=1 THEORY WILL BE NORMALIZED AND RETURN IS MADE WITH IT STILL
C    NORMALIZED
C    KND=2 THEORY WILL BE NORMALIZED BUT RETURN IS MADE WITH IT RESTORED
C
C---- FIT AVN BY A1+A2*X+EXP(A3*X)*(A4*COST**3 - A5*COST**2 + A6*COST)
C
C---- SEARCH THE MIN OF THE AVN
      AVMIN = 100000.
        DO 4 K=1,NPTS
    4   IF(AVN(K).LT.AVMIN) AVMIN = AVN(K)
C
C---- PRESET PARAMETERS
      A(1) = AVMIN-1.0
      A(2) = 1.5*COST+0.5
      A(3) = AEXPA
      A(4) = (1. + COST) * .01
C
C     DO 980  J = 1,3
C     DELTAA(J) = .2 * A(J)
C80   IF(DELTAA(J).GT.5.) DELTAA(J) = 5.
C
      FLAMDA = 0.001D0
      CHISQO = 0.
      NPASS = 0
   20 CONTINUE
      CALL CURFIT(X,AVN,ERR,NPTS,NTERMS,MODE,
     1          A,DELTAA,SIGMAA,FLAMDA,YFIT,CHISQ,IPRO)
C    1          A,DELTAA,SIGMAA,FLAMDA,YFIT,CHISQ)
C---- CHECK CONVERGENCE
      DCHISQ = CHISQ-CHISQO
      DCHISQ = ABS(DCHISQ)/CHISQ
      CHISQO = CHISQ
C
C******** TEST DUMP
      IF(IPRO.EQ.1)  WRITE(6,6666) CHISQO,DCHISQ,(A(K),K=1,NTERMS)
 6666 FORMAT(' CHI,D: ',1P2E11.3,' A: ',8E11.3)
C
      NPASS = NPASS+1
      IF(NPASS.GT.100) GO TO 90
      IF(DCHISQ.GT.0.0005) GO TO 20
   90 CONTINUE
C
C     CALL HISTOG(20,SD(1,ND),YFIT,0.,1.,2,4,DUMMY)
C
      DO 1232  JJ = 1,NPTS
      WGT = SD(JJ,ND)
      CALL HFILL(ND,X(JJ),0,WGT)
      WGT = YFIT(JJ)
      CALL HFILL(30+ND,X(JJ),0,WGT)
1232  CONTINUE
      WGT = CC1
      CALL HFILL(151,COST,0,WGT)
      WGT = CC
      CALL HFILL(152,COST,0,WGT)
      WGT = PART1
      CALL HFILL(153,COST,0,WGT)
      WGT = PART2
      CALL HFILL(154,COST,0,WGT)
      WGT = PART3
      CALL HFILL(155,COST,0,WGT)
C
C---- PRINT INPUT
C
      WRITE(6,600) COST,AVN,ERR
  600 FORMAT('0 COST(ORIGINE)=',F8.2,/
     1' AVN ',3X,2(5F7.2,2X),/,' ',2(5F7.2,2X),
     2      /' +/-',2(5F7.2,2X),/,' ',2(5F7.2,2X))
C
      WRITE(6,602) NPASS
  602 FORMAT(' TRIAL NO.(MAX=100)=',I4)
C
   22 WRITE(6,601) (YFIT(L),L=1,NPTS),CHISQO,CHISQ,DCHISQ,
     1      (A(K),SIGMAA(K),K=1,NTERMS)
  601 FORMAT(' YFIT:',2(5F7.1,2X),/,' ',2(5F7.1,2X),/' CHISQO,...=',
     1 3E12.4,/,' A:',1P4(E11.3,'(',E11.3,')',2X))
   23 IF(ND.GT.50) GO TO 1
      RESLT(ND,1) = COST
      RESLT(ND,2) = A(1)
      RESLT(ND,3) = SIGMAA(1)
      RESLT(ND,4) = A(2)
      RESLT(ND,5) = SIGMAA(2)
      RESLT(ND,6) = A(3)
      RESLT(ND,7) = SIGMAA(3)
      RESLT(ND,8) = A(4)
      RESLT(ND,9) = SIGMAA(4)
      GO TO 1
  100 WRITE(6,603)
  603 FORMAT(' *************** FIT RESULTS **********')
        DO 30 N=1,ND
        WRITE(6,604) (RESLT(N,K),K=1,9)
  604   FORMAT(' FITRES ',F6.2,4(F11.3,' +/-',F6.3))
   30   CONTINUE
C
C******************************************************
C
C    FIT A(1) - A(3) AS FUNCTIONS OF COST
C
      IPRO = 0
C  STEP IN COST
      DCST = ABS(RESLT(1,1)-RESLT(2,1))
C
      DO 300 MODEFT=2,5
      NTERMS = 8
C
C   THERE ARE ND DIFFERENT COST POINTS
C
      NPTS = ND
C
        DO 400 K=1,8
  400   A(K) = 0.
      IF(MODEFT.NE.2) GO TO 403
      A(1) = 22.6
      A(2) = -1.9
      A(3) = -285.0
      A(4) = 1459.
      A(5) = -2780.0
      A(6) = 2306.0
      A(7) = -690.0
      GO TO 404
  403 CONTINUE
      IF(MODEFT.NE.3) GO TO 433
      A(1) = 19.95
      A(2) = 1.58
      A(3) = -0.3
      A(4) = -2.0
      A(6) = 1.45
      A(7) = 0.5
      GO TO 404
  433 CONTINUE
      IF(MODEFT.NE.4) GO TO 434
      A(1) = AEXPA
      GO TO 404
  434 CONTINUE
      A(1) = 0.01
      A(2) = 0.01
  404 CONTINUE
      FLAMDA = 0.001D0
      CHISQO = 0.
      NPASS = 0
      IPOS = (MODEFT-1)*2
      IPOS1 = IPOS+1
  200 CONTINUE
      IF(NPASS.EQ.1) WRITE(6,1471) MODEFT
1471  FORMAT(' FIT2:  MODEFT = ',I3)
      CALL CURFIT(RESLT(1,1),RESLT(1,IPOS),RESLT(1,IPOS1),
     1      NPTS,NTERMS,MODE,A,DELTAA,SIGMAA,FLAMDA,
     2      RESLT(1,MODEFT+8),CHISQ,IPRO)
C    2      RESLT(1,MODEFT+8),CHISQ)
C---- CHECK CONVERGENCE
C     IF(NPASS.EQ.1)
C    1CALL HISTOG(20,RESLT(1,IPOS),RESLT(1,MODEFT+8),1.,-DCST,2,4,DUMMY)
      WRITE(6,1781) NPASS,CHISQ
1781  FORMAT(' FIT2: NPASS CHISQ ',I4,E14.6)
      DCHISQ = CHISQ-CHISQO
      DCHISQ = ABS(DCHISQ)/CHISQ
C
      WRITE(6,6666) CHISQ,DCHISQ,(A(K),K=1,NTERMS)
      CHISQO = CHISQ
      NPASS = NPASS+1
      IF(NPASS.GT.100) GO TO 250
      IF(DCHISQ.GT.0.0005) GO TO 200
C---- PRINT FIT PARAMETERS
C
C GRAPHICAL DISPLAY OF INPUT AND FIT
C
C     CALL HISTOG(20,RESLT(1,IPOS),RESLT(1,MODEFT+8),1.,-DCST,2,4,DUMMY)
C
      COSS = 1.+DCST
      DO 2132  JJ = 1,ND
      COSS = COSS - DCST
      WGT = RESLT(JJ,IPOS)
      CALL HFILL(99+MODEFT,COSS,0,WGT)
      WGT = RESLT(JJ,MODEFT+8)
      CALL HFILL(109+MODEFT,COSS,0,WGT)
2132  CONTINUE
C
  250 WRITE(6,602) NPASS
      WRITE(6,605) A,SIGMAA
  605 FORMAT('0  **** A ***',1P8E12.4,/'       +/-   ',8E12.4)
C
        DO 260 K=1,8
  260   BC(K,MODEFT-1) = A(K)
  300 CONTINUE
C
        DO 220 N=1,ND
        WRITE(6,606) (RESLT(N,K),K=1,13)
  606   FORMAT(' RESLT:',F8.2,12F9.3)
  220   CONTINUE
C
C ************ ALL FITS ARE FINISHED *********
C
C    PRINT THE FINAL RESULT
C
       WRITE(6,2271)
2271   FORMAT(1H1)
C
        DO 1000 N = 1,ND
        CHISQ = 0.
C-----   SET PARAMETERS
          DO 1050 IP=1,4
          MODEFT = IP+1
          A(IP) = FUNCTN(RESLT(1,1),N,BC(1,IP))
 1050     CONTINUE
        MODEFT = 1
        COST = RESLT(N,1)
        IXE = 20
        IF(IENDCP.EQ.1) IXE = 12
          DO 1100 IX=1,IXE
          YFIT(IX) = FUNCTN(X,IX,A)
          CHISQ = CHISQ+((SD(IX,N)-YFIT(IX))/SER(IX,N))**2
 1100     CONTINUE
C
          DO 2491  IX = 1,IXE
          WGT = YFIT(IX)
          CALL HFILL(60+N,X(IX),0,WGT)
2491  CONTINUE
          CALL HFILL(121,COST,0,CHISQ)
          CALL HFILL(161,COST,0,CC1)
          CALL HFILL(162,COST,0,CC)
          CALL HFILL(163,COST,0,PART1)
          CALL HFILL(164,COST,0,PART2)
          CALL HFILL(165,COST,0,PART3)
C
        WRITE(6,607) COST,CHISQ,A(1),A(2),A(3),A(4),
     1           (SD(K,N),K=1,IXE),(YFIT(L),L=1,IXE)
  607   FORMAT('0COST=',F6.3,3X,'CHISQ=',E10.3,' A1-A4=',
     1        4F10.4,/' SD ',10F7.2,/,' ',10F7.2,
     2        /' YFIT ',10F7.2,/,' ',10F7.2)
 1000   CONTINUE
      CALL HISTDO
      STOP
      END
C***** CERENKOV LIGHT DETECTION EFFICIENCY FIT *****
C     S.YAMADA   15-12-77
C     RE-CODED ON 15-12-77 TO RECOVER LOST VERSION
C     UPDATED  2-10-85 TO SF6, J.OLSSON
C     FORMATS MODIFIED, 5.1.1986    J.OLSSON
C  2 COMPONENT VERSION 30.1.1986    J.OLSSON
C  A4-A7 IN CC READ IN, VARIABLE
C
C         GEP VERSION
C
C     READ CONVERTED PHOTOELECTRON NUMBERS FROM THE TABLE
C     FOR DIFFERENT DIRECTION AND POSITIONS OF THE ORIGINS.
C---- FIT THE NUMBERS FOR FIXED DIRECTION FIRST AND THEN THE FIT
C     PARAMETERS ARE FITTED AS A FUNCTION OF THE DIRECTION.
C
C   SEE DESCRIPTION IN JADE NOTE 20, BY S.YAMADA
C
C
      REAL*8 FLAMDA
      DIMENSION AVN(20),ERR(20), TRY(20),RESLT(50,13)
      DIMENSION X(20),YFIT(21), A(8),DELTAA(8),SIGMAA(8)
      DATA AVN/20*0./,ERR/20*0./,TRY/19*100.,60./,DELTAA/8*1.0/
      DATA X/0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5,
     1  13.5,14.5,15.5,16.5,17.5,18.5,19.3/
      DATA MODE/1/, ND/0/
      DIMENSION SD(20,50), SER(20,50), BC(8,4)
      LOGICAL*1 TEXT
      DIMENSION TEXT(72)
C
      COMMON /COSTB/ COST,MODEFT,AESS,BASI,CC,CC1,AEXPA,A4,A5,A6,A7
     $ ,PART1,PART2,PART3
C
      DATA IENDCP/0/
C
C---- IENDCP=1  FOR THE ENDCAP, IENDCP=0  FOR THE BARREL PART.
C
      MODEFT = 1
      NTERMS = 4
      IPRO = 0
      NPTS = 20
      IF(IENDCP.EQ.1) NPTS = 12
      IF(IENDCP.EQ.1) X(12) = 11.1
C
C PLOT INITIALIZATION
C
      CALL ACORE(10000)
      DO 3705  I = 1,21
      CALL HBOOK1( I,' NR OF PHOTONS VS. RAD. LENGTH$',20,0.,20.)
      CALL HBOOK1(30+ I,' NR OF PHOTONS FITTED, SINGLE$',20,0.,20.)
      CALL HBOOK1(60+ I,' NR OF PHOTONS FITTED, COMBINED$',20,0.,20.)
3705  CONTINUE
      DO 3706  I = 1,4
      CALL HBOOK1(100+ I,' FIT PARAMETER VS. COSTHETA$',21,-1.05,1.05)
      CALL HBOOK1(110+ I,' FIT PARAMETER, GLOBAL FIT$',21,-1.05,1.05)
3706  CONTINUE
      CALL HBOOK1(121,'CHISQ VS. COSTHETA$',21,-1.05,1.05)
      CALL HBOOK1(151,' CC1(COST)VS COSTHETA$',21,-1.05,1.05)
      CALL HBOOK1(152,' CC=10**CC1 VS COSTHETA$',21,-1.05,1.05)
      CALL HBOOK1(153,' PART1 VS COSTHETA$',21,-1.05,1.05)
      CALL HBOOK1(154,' PART2 VS COSTHETA$',21,-1.05,1.05)
      CALL HBOOK1(155,' PART3 VS COSTHETA$',21,-1.05,1.05)
      CALL HBOOK1(161,' CC1(COST)VS COSTHETA,GLOBAL$',21,-1.05,1.05)
      CALL HBOOK1(162,' CC=10**CC1 VS COSTHETA,GLOBAL$',21,-1.05,1.05)
      CALL HBOOK1(163,' PART1 VS COSTHETA,GLOBAL$',21,-1.05,1.05)
      CALL HBOOK1(164,' PART2 VS COSTHETA,GLOBAL$',21,-1.05,1.05)
      CALL HBOOK1(165,' PART3 VS COSTHETA,GLOBAL$',21,-1.05,1.05)
C
C
7608  FORMAT(72A1)
7508  FORMAT(' ',72A1)
      READ 7608, TEXT
      WRITE(6,7508) TEXT
C
C---- READ THE MONTE CARLO RESULTS.
C---- AESS,BASI,AEXPA,A4-7 ARE PARAMETERS USED IN FUNCTN
C
      READ(5,501) AESS,BASI,AEXPA,A4,A5,A6,A7
  501 FORMAT(3F10.2,4F10.4)
      WRITE(6,6610) AESS,BASI,AEXPA,A4,A5,A6,A7
6610  FORMAT('  FUNCTN PAR.: AESS BASI AEXPA ',3F10.2,' A4-7 ',4F10.4)
C
C  LOOP OVER VARIOUS COST STARTING HERE   ###########################
C
 1    CONTINUE
      READ(5,500,END=100) COST
      READ(5,502,END=100) (AVN(K),K=1,10)
      READ(5,502,END=100) (AVN(K),K=11,NPTS)
  500 FORMAT(F6.3)
  502 FORMAT(10F7.2)
C  COMPENSATE FOR INCOMPLETE STEP AT END
C     AVN(NPTS) = AVN(NPTS)*1.7
C
C---- CALCULATE THE STATISTICAL ERRORS
C---- FORMULA ASSUMES AVN IS MEAN OF TRY ATTEMPTS, ERR IS THEN
C---- ERROR ON MEAN:  SIGMA**2/TRY, WITH SIGMA = ROOT(AVN)
C
        DO 2 K=1,NPTS
    2   ERR(K) = SQRT(AVN(K)/TRY(K))*0.5
C----   THE FACTOR 0.5 IS DUE TO THE CONVERSION EFFICIENCY MAX.(0.25).
C
C  SINCE AVN IS NOT A MEAN, BUT THE SUM OF TRY ATTEMPTS, THE CORRECTER
C  ERROR IS SIMPLY SQRT(SD)
C
C       DO 2 K=1,NPTS
C   2   ERR(K) = SQRT(AVN(K))
C
C---- STORE THE SOURCE DATA FOR FUTURE USE.
      ND = ND+1
      IF(ND.GT.50) GO TO 33
        DO 3 K=1,NPTS
        SD(K,ND) = AVN(K)
    3   SER(K,ND) = ERR(K)
33     CONTINUE
C
C
C PLOT INPUT
C
C     CALL HISTOG(20,SD(1,ND),X,0.,1.,0,4,DUMMY)
C     SUBROUTINE HISTOG (NBIN,CINCON,THEORY,BINLOW,BINSIZ,KND,NL,HISTIT)
C    NBIN = NO. OF BINS
C    CINCON = EXPERIMENTAL DISTRIBUTION
C    THEORY=THEORETICAL DISTRIBUTION (SUPERIMPOSED IF REQUIRED)
C    BINLOW = LOW EDGE OF FIRST BIN
C    BINSIZ = SIZE OF BIN
C    KND=0 IF YOU DON'T WISH TO SUPPLY A THEORETICAL CURVE
C    KND=1 THEORY WILL BE NORMALIZED AND RETURN IS MADE WITH IT STILL
C    NORMALIZED
C    KND=2 THEORY WILL BE NORMALIZED BUT RETURN IS MADE WITH IT RESTORED
C
C---- FIT AVN BY A1+A2*X+EXP(A3*X)*(A4*COST**3 - A5*COST**2 + A6*COST)
C
C---- SEARCH THE MIN OF THE AVN
      AVMIN = 100000.
        DO 4 K=1,NPTS
    4   IF(AVN(K).LT.AVMIN) AVMIN = AVN(K)
C
C---- PRESET PARAMETERS
      A(1) = AVMIN-1.0
      A(2) = 1.5*COST+0.5
      A(3) = AEXPA
      A(4) = (1. + COST) * .01
C
C     DO 980  J = 1,3
C     DELTAA(J) = .2 * A(J)
C80   IF(DELTAA(J).GT.5.) DELTAA(J) = 5.
C
      FLAMDA = 0.001D0
      CHISQO = 0.
      NPASS = 0
   20 CONTINUE
      CALL CURFIT(X,AVN,ERR,NPTS,NTERMS,MODE,
     1          A,DELTAA,SIGMAA,FLAMDA,YFIT,CHISQ,IPRO)
C    1          A,DELTAA,SIGMAA,FLAMDA,YFIT,CHISQ)
C---- CHECK CONVERGENCE
      DCHISQ = CHISQ-CHISQO
      DCHISQ = ABS(DCHISQ)/CHISQ
      CHISQO = CHISQ
C
C******** TEST DUMP
      IF(IPRO.EQ.1)  WRITE(6,6666) CHISQO,DCHISQ,(A(K),K=1,NTERMS)
 6666 FORMAT(' CHI,D: ',1P2E11.3,' A: ',8E11.3)
C
      NPASS = NPASS+1
      IF(NPASS.GT.100) GO TO 90
      IF(DCHISQ.GT.0.0005) GO TO 20
   90 CONTINUE
C
C     CALL HISTOG(20,SD(1,ND),YFIT,0.,1.,2,4,DUMMY)
C
      DO 1232  JJ = 1,NPTS
      WGT = SD(JJ,ND)
      CALL HFILL(ND,X(JJ),0,WGT)
      WGT = YFIT(JJ)
      CALL HFILL(30+ND,X(JJ),0,WGT)
1232  CONTINUE
      WGT = CC1
      CALL HFILL(151,COST,0,WGT)
      WGT = CC
      CALL HFILL(152,COST,0,WGT)
      WGT = PART1
      CALL HFILL(153,COST,0,WGT)
      WGT = PART2
      CALL HFILL(154,COST,0,WGT)
      WGT = PART3
      CALL HFILL(155,COST,0,WGT)
C
C---- PRINT INPUT
C
      WRITE(6,600) COST,AVN,ERR
  600 FORMAT('0 COST(ORIGINE)=',F8.2,/
     1' AVN ',3X,2(5F7.2,2X),/,' ',2(5F7.2,2X),
     2      /' +/-',2(5F7.2,2X),/,' ',2(5F7.2,2X))
C
      WRITE(6,602) NPASS
  602 FORMAT(' TRIAL NO.(MAX=100)=',I4)
C
   22 WRITE(6,601) (YFIT(L),L=1,NPTS),CHISQO,CHISQ,DCHISQ,
     1      (A(K),SIGMAA(K),K=1,NTERMS)
  601 FORMAT(' YFIT:',2(5F7.1,2X),/,' ',2(5F7.1,2X),/' CHISQO,...=',
     1 3E12.4,/,' A:',1P4(E11.3,'(',E11.3,')',2X))
   23 IF(ND.GT.50) GO TO 1
      RESLT(ND,1) = COST
      RESLT(ND,2) = A(1)
      RESLT(ND,3) = SIGMAA(1)
      RESLT(ND,4) = A(2)
      RESLT(ND,5) = SIGMAA(2)
      RESLT(ND,6) = A(3)
      RESLT(ND,7) = SIGMAA(3)
      RESLT(ND,8) = A(4)
      RESLT(ND,9) = SIGMAA(4)
      GO TO 1
  100 WRITE(6,603)
  603 FORMAT(' *************** FIT RESULTS **********')
        DO 30 N=1,ND
        WRITE(6,604) (RESLT(N,K),K=1,9)
  604   FORMAT(' FITRES ',F6.2,4(F11.3,' +/-',F6.3))
   30   CONTINUE
C
C******************************************************
C
C    FIT A(1) - A(3) AS FUNCTIONS OF COST
C
      IPRO = 0
C  STEP IN COST
      DCST = ABS(RESLT(1,1)-RESLT(2,1))
C
      DO 300 MODEFT=2,5
      NTERMS = 8
C
C   THERE ARE ND DIFFERENT COST POINTS
C
      NPTS = ND
C
        DO 400 K=1,8
  400   A(K) = 0.
      IF(MODEFT.NE.2) GO TO 403
      A(1) = 22.6
      A(2) = -1.9
      A(3) = -285.0
      A(4) = 1459.
      A(5) = -2780.0
      A(6) = 2306.0
      A(7) = -690.0
      GO TO 404
  403 CONTINUE
      IF(MODEFT.NE.3) GO TO 433
      A(1) = 19.95
      A(2) = 1.58
      A(3) = -0.3
      A(4) = -2.0
      A(6) = 1.45
      A(7) = 0.5
      GO TO 404
  433 CONTINUE
      IF(MODEFT.NE.4) GO TO 434
      A(1) = AEXPA
      GO TO 404
  434 CONTINUE
      A(1) = 0.01
      A(2) = 0.01
  404 CONTINUE
      FLAMDA = 0.001D0
      CHISQO = 0.
      NPASS = 0
      IPOS = (MODEFT-1)*2
      IPOS1 = IPOS+1
  200 CONTINUE
      IF(NPASS.EQ.1) WRITE(6,1471) MODEFT
1471  FORMAT(' FIT2:  MODEFT = ',I3)
      CALL CURFIT(RESLT(1,1),RESLT(1,IPOS),RESLT(1,IPOS1),
     1      NPTS,NTERMS,MODE,A,DELTAA,SIGMAA,FLAMDA,
     2      RESLT(1,MODEFT+8),CHISQ,IPRO)
C    2      RESLT(1,MODEFT+8),CHISQ)
C---- CHECK CONVERGENCE
C     IF(NPASS.EQ.1)
C    1CALL HISTOG(20,RESLT(1,IPOS),RESLT(1,MODEFT+8),1.,-DCST,2,4,DUMMY)
      WRITE(6,1781) NPASS,CHISQ
1781  FORMAT(' FIT2: NPASS CHISQ ',I4,E14.6)
      DCHISQ = CHISQ-CHISQO
      DCHISQ = ABS(DCHISQ)/CHISQ
C
      WRITE(6,6666) CHISQ,DCHISQ,(A(K),K=1,NTERMS)
      CHISQO = CHISQ
      NPASS = NPASS+1
      IF(NPASS.GT.100) GO TO 250
      IF(DCHISQ.GT.0.0005) GO TO 200
C---- PRINT FIT PARAMETERS
C
C GRAPHICAL DISPLAY OF INPUT AND FIT
C
C     CALL HISTOG(20,RESLT(1,IPOS),RESLT(1,MODEFT+8),1.,-DCST,2,4,DUMMY)
C
      COSS = 1.+DCST
      DO 2132  JJ = 1,ND
      COSS = COSS - DCST
      WGT = RESLT(JJ,IPOS)
      CALL HFILL(99+MODEFT,COSS,0,WGT)
      WGT = RESLT(JJ,MODEFT+8)
      CALL HFILL(109+MODEFT,COSS,0,WGT)
2132  CONTINUE
C
  250 WRITE(6,602) NPASS
      WRITE(6,605) A,SIGMAA
  605 FORMAT('0  **** A ***',1P8E12.4,/'       +/-   ',8E12.4)
C
        DO 260 K=1,8
  260   BC(K,MODEFT-1) = A(K)
  300 CONTINUE
C
        DO 220 N=1,ND
        WRITE(6,606) (RESLT(N,K),K=1,13)
  606   FORMAT(' RESLT:',F8.2,12F9.3)
  220   CONTINUE
C
C ************ ALL FITS ARE FINISHED *********
C
C    PRINT THE FINAL RESULT
C
       WRITE(6,2271)
2271   FORMAT(1H1)
C
        DO 1000 N = 1,ND
        CHISQ = 0.
C-----   SET PARAMETERS
          DO 1050 IP=1,4
          MODEFT = IP+1
          A(IP) = FUNCTN(RESLT(1,1),N,BC(1,IP))
 1050     CONTINUE
        MODEFT = 1
        COST = RESLT(N,1)
        IXE = 20
        IF(IENDCP.EQ.1) IXE = 12
          DO 1100 IX=1,IXE
          YFIT(IX) = FUNCTN(X,IX,A)
          CHISQ = CHISQ+((SD(IX,N)-YFIT(IX))/SER(IX,N))**2
 1100     CONTINUE
C
          DO 2491  IX = 1,IXE
          WGT = YFIT(IX)
          CALL HFILL(60+N,X(IX),0,WGT)
2491  CONTINUE
          CALL HFILL(121,COST,0,CHISQ)
          CALL HFILL(161,COST,0,CC1)
          CALL HFILL(162,COST,0,CC)
          CALL HFILL(163,COST,0,PART1)
          CALL HFILL(164,COST,0,PART2)
          CALL HFILL(165,COST,0,PART3)
C
        WRITE(6,607) COST,CHISQ,A(1),A(2),A(3),A(4),
     1           (SD(K,N),K=1,IXE),(YFIT(L),L=1,IXE)
  607   FORMAT('0COST=',F6.3,3X,'CHISQ=',E10.3,' A1-A4=',
     1        4F10.4,/' SD ',10F7.2,/,' ',10F7.2,
     2        /' YFIT ',10F7.2,/,' ',10F7.2)
 1000   CONTINUE
      CALL HISTDO
      STOP
      END
C   26/06/78 601301806  MEMBER NAME  @CERSF55 (S)           FORTRAN
C***** CERENKOV LIGHT DETECTION EFFICIENCY FIT *****
C     S.YAMADA   15-12-77
C     RE-CODED ON 15-12-77 TO RECOVER LOST VERSION
C     UPDATED  2-10-85 TO SF6, J.OLSSON
C     FORMATS MODIFIED, 5.1.1986    J.OLSSON
C   VERSION TO FIT SF5 AS TEST   7.1.86  J.OLSSON
C   MODIFIED BACK TO @CERFIT1 VERSION ON F22YAM.WORKS     29.1.1986
C
C     READ CONVERTED PHOTOELECTRON NUMBERS FROM THE TABLE
C     FOR DIFFERENT DIRECTION AND POSITIONS OF THE ORIGINS.
C---- FIT THE NUMBERS FOR FIXED DIRECTION FIRST AND THEN THE FIT
C     PARAMETERS ARE FITTED AS A FUNCTION OF THE DIRECTION.
C
C   SEE DESCRIPTION IN JADE NOTE 20, BY S.YAMADA
C
C
      REAL*8 FLAMDA
      DIMENSION AVN(14),ERR(14), TRY(14),RESLT(50,8)
      DIMENSION X(14),YFIT(14), A(8),DELTAA(8),SIGMAA(8)
      DATA AVN/14*0./,ERR/14*0./,TRY/13*20.,8./,DELTAA/8*1.0/
      DATA X/0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5,
     1  13.2/
      DATA MODE/1/, ND/0/
      DIMENSION SD(14,50), SER(14,50), BC(8,2)
      LOGICAL*1 TEXT
      DIMENSION TEXT(72)
C
      COMMON /COSTB/ COST,MODEFT
C
      DATA IENDCP/0/
C
C---- IENDCP=1  FOR THE ENDCAP, IENDCP=0  FOR THE BARREL PART.
C
      MODEFT = 1
      NTERMS = 2
      IPRO = 0
      NPTS = 14
      IF(IENDCP.EQ.1) NPTS = 12
      IF(IENDCP.EQ.1) X(12) = 11.1
C
C PLOT INITIALIZATION
C
      CALL ACORE(10000)
      DO 3705  I = 1,21
      CALL HBOOK1( I,' NR OF PHOTONS VS. RAD. LENGTH$',20,0.,20.)
      CALL HBOOK1(30+ I,' NR OF PHOTONS FITTED, SINGLE$',20,0.,20.)
      CALL HBOOK1(60+ I,' NR OF PHOTONS FITTED, COMBINED$',20,0.,20.)
3705  CONTINUE
      DO 3706  I = 1,2
      CALL HBOOK1(100+ I,' FIT PARAMETER VS. COSTHETA$',21,-1.05,1.05)
      CALL HBOOK1(110+ I,' FIT PARAMETER, GLOBAL FIT$',21,-1.05,1.05)
3706  CONTINUE
      CALL HBOOK1(121,'CHISQ VS. COSTHETA$',21,-1.05,1.05)
C
C
7608  FORMAT(72A1)
7508  FORMAT(' ',72A1)
      READ 7608, TEXT
      WRITE(6,7508) TEXT
C
C---- READ THE MONTE CARLO RESULTS.
C---- IPR:  RAW DATA DUMP
C---- IPRR:  DUMP OF FIT RESULTS
C
      READ(5,501) IPR,IPRR
  501 FORMAT(2I10)
C
C  LOOP OVER VARIOUS COST STARTING HERE   ###########################
C
 1    CONTINUE
      READ(5,500,END=100) COST
      READ(5,502,END=100) (AVN(K),K=1,10)
      READ(5,502,END=100) (AVN(K),K=11,NPTS)
  500 FORMAT(F6.3)
  502 FORMAT(10F7.2)
C  COMPENSATE FOR INCOMPLETE STEP AT END
C    13.2 IS CENTER, ==> 40 %  ==>  FACTOR 2.5 TO COMPENSATE
C     AVN(NPTS) = AVN(NPTS)*2.5
C
C---- CALCULATE THE STATISTICAL ERRORS
C---- FORMULA ASSUMES AVN IS MEAN OF TRY ATTEMPTS, ERR IS THEN
C---- ERROR ON MEAN:  SIGMA**2/TRY
C
        DO 2 K=1,NPTS
    2   ERR(K) = SQRT(AVN(K)/TRY(K))*0.5
C----   THE FACTOR 0.5 IS DUE TO THE CONVERSION EFFICIENCY MAX.(0.25).
C
C  SINCE AVN IS NOT A MEAN, BUT THE SUM OF TRY ATTEMPTS, THE CORRECTER
C  ERROR IS SIMPLY SQRT(SD)
C
C       DO 2 K=1,NPTS
C   2   ERR(K) = SQRT(AVN(K))
C
C---- STORE THE SOURCE DATA FOR FUTURE USE.
      ND = ND+1
      IF(ND.GT.50) GO TO 33
        DO 3 K=1,NPTS
        SD(K,ND) = AVN(K)
    3   SER(K,ND) = ERR(K)
33     CONTINUE
C
C
C PLOT INPUT
C
C     CALL HISTOG(20,SD(1,ND),X,0.,1.,0,4,DUMMY)
C     SUBROUTINE HISTOG (NBIN,CINCON,THEORY,BINLOW,BINSIZ,KND,NL,HISTIT)
C    NBIN = NO. OF BINS
C    CINCON = EXPERIMENTAL DISTRIBUTION
C    THEORY=THEORETICAL DISTRIBUTION (SUPERIMPOSED IF REQUIRED)
C    BINLOW = LOW EDGE OF FIRST BIN
C    BINSIZ = SIZE OF BIN
C    KND=0 IF YOU DON'T WISH TO SUPPLY A THEORETICAL CURVE
C    KND=1 THEORY WILL BE NORMALIZED AND RETURN IS MADE WITH IT STILL
C    NORMALIZED
C    KND=2 THEORY WILL BE NORMALIZED BUT RETURN IS MADE WITH IT RESTORED
C
C---- FIT AVN BY A1+A2*X+EXP(A3*X)*(A4*COST**3 - A5*COST**2 + A6*COST)
C
C---- SEARCH THE MIN OF THE AVN
      AVMIN = 1000.
        DO 4 K=1,NPTS
    4   IF(AVN(K).LT.AVMIN) AVMIN = AVN(K)
C
C---- PRESET PARAMETERS
      A(1) = AVMIN-1.0
      A(2) = 1.5*COST+0.5
      A(3) = 0.001
      A(4) = 0.55
C
C     DO 980  J = 1,2
C     DELTAA(J) = .2 * A(J)
C80   IF(DELTAA(J).GT.5.) DELTAA(J) = 5
C
      FLAMDA = 0.001D0
C     FLAMDA = 0.001
      CHISQO = 0.
      NPASS = 0
   20 CONTINUE
      CALL CURFIT(X,AVN,ERR,NPTS,NTERMS,MODE,
     1          A,DELTAA,SIGMAA,FLAMDA,YFIT,CHISQ,IPRO)
C    1          A,DELTAA,SIGMAA,FLAMDA,YFIT,CHISQ)
C---- CHECK CONVERGENCE
      DCHISQ = CHISQ-CHISQO
      DCHISQ = ABS(DCHISQ)/CHISQ
      CHISQO = CHISQ
C
C******** TEST DUMP
      IF(IPRO.EQ.1)  WRITE(6,6666) CHISQO,DCHISQ,(A(K),K=1,NTERMS)
 6666 FORMAT(' CHI,D: ',1P2E11.3,' A: ',8E11.3)
C
      NPASS = NPASS+1
      WRITE(6,602) NPASS
      IF(NPASS.GT.100) GO TO 90
      IF(DCHISQ.GT.0.0005) GO TO 20
   90 CONTINUE
C
C     CALL HISTOG(14,SD(1,ND),YFIT,0.,1.,2,4,DUMMY)
C
      DO 1232  JJ = 1,NPTS
      WGT = SD(JJ,ND)
      CALL HFILL(ND,X(JJ),0,WGT)
      WGT = YFIT(JJ)
      CALL HFILL(30+ND,X(JJ),0,WGT)
1232  CONTINUE
C---- PRINT INPUT
C
      IF(IPR.NE.0) WRITE(6,600) COST,AVN,ERR
  600 FORMAT('0 COST(ORIGINE)=',F8.2,/
     1' AVN ',3X,2(5F7.2,2X),/,' ',4F7.2,
     2      /' +/-',2(5F7.2,2X),/,' ',4F7.2)
C
      WRITE(6,602) NPASS
  602 FORMAT(' TRIAL NO.(MAX=100)=',I4)
C
   22 WRITE(6,601) YFIT,CHISQO,CHISQ,DCHISQ,
     1      (A(K),SIGMAA(K),K=1,NTERMS)
  601 FORMAT(' YFIT:',2(5F7.1,2X),/,' ',4F7.1,/' CHISQO,...=',
     1 3E12.4,/,' A:',1P4(E11.3,'(',E11.3,')',2X))
   23 IF(ND.GT.50) GO TO 1
      RESLT(ND,1) = COST
      RESLT(ND,2) = A(1)
      RESLT(ND,3) = SIGMAA(1)
      RESLT(ND,4) = A(2)
      RESLT(ND,5) = SIGMAA(2)
      RESLT(ND,6) = A(3)
      RESLT(ND,7) = SIGMAA(3)
      RESLT(ND,8) = A(4)
      GO TO 1
  100 WRITE(6,603)
  603 FORMAT(' *************** FIT RESULTS **********')
        DO 30 N=1,ND
        WRITE(6,604) (RESLT(N,K),K=1,8)
  604   FORMAT(' FITRES ',F6.2,3(F11.3,' +/-',F6.3),F7.3)
   30   CONTINUE
C
C******************************************************
C
C    FIT A(1) - A(2) AS FUNCTIONS OF COST
C
      IPRO = 0
C  STEP IN COST
      DCST = ABS(RESLT(1,1)-RESLT(2,1))
C
      DO 300 MODEFT=2,3
      NTERMS = 8
C
C   THERE ARE ND DIFFERENT COST POINTS
C
      NPTS = ND
C
        DO 400 K=1,8
  400   A(K) = 0.
      IF(MODEFT.EQ.3) GO TO 403
C     A(1) = 22.6
C     A(2) = -1.9
C     A(3) = -285.0
C     A(4) = 1459.
C     A(5) = -2780.0
C     A(6) = 2306.0
C     A(7) = -690.0
      A(1) = 58.22
      A(2) = 29.2
      A(3) = -1150.6
      A(4) = 4664.1
      A(5) = -6253.8
      A(6) = 685.51
      A(7) = 4351.7
      A(8) = -2296.5
      GO TO 404
 403  CONTINUE
C     A(1) = 19.95
C     A(2) = 1.58
C     A(3) = -0.3
C     A(4) = -2.0
C     A(6) = 1.45
C     A(7) = 0.5
      A(1) = 0.4225
      A(2) = 1.0556
      A(3) = -1.6712
      A(4) = 9.0240
      A(5) = 2.8498
      A(6) = -23.35
      A(7) = -0.9666
      A(8) = 15.58
  404 CONTINUE
      FLAMDA = 0.001D0
C     FLAMDA = 0.001
      CHISQO = 0.
      NPASS = 0
      IPOS = (MODEFT-1)*2
      IPOS1 = IPOS+1
  200 CONTINUE
      CALL CURFIT(RESLT(1,1),RESLT(1,IPOS),RESLT(1,IPOS1),
     1      NPTS,NTERMS,MODE,A,DELTAA,SIGMAA,FLAMDA,
     2      RESLT(1,MODEFT+4),CHISQ,IPRO)
C    2      RESLT(1,MODEFT+4),CHISQ)
C---- CHECK CONVERGENCE
C     IF(NPASS.EQ.1) CALL
C    1 HISTOG(ND,RESLT(1,IPOS),RESLT(1,MODEFT+6),1.05,-DCST,2,4,DUMMY)
      IF(NPASS.EQ.1) WRITE(6,1471) MODEFT
1471  FORMAT(' FIT2:  MODEFT = ',I3)
      WRITE(6,1781) NPASS,CHISQ
1781  FORMAT(' FIT2: NPASS CHISQ ',I4,E14.6)
      DCHISQ = CHISQ-CHISQO
      DCHISQ = ABS(DCHISQ)/CHISQ
C
      WRITE(6,6666) CHISQ,DCHISQ,(A(K),K=1,NTERMS)
      CHISQO = CHISQ
      NPASS = NPASS+1
      IF(NPASS.GT.100) GO TO 250
      IF(DCHISQ.GT.0.0005) GO TO 200
C---- PRINT FIT PARAMETERS
C
C GRAPHICAL DISPLAY OF INPUT AND FIT
C
C     CALL
C    $ HISTOG(ND,RESLT(1,IPOS),RESLT(1,MODEFT+4),1.05,-DCST,2,4,DUMMY)
C
      COSS = 1.+DCST
      DO 2132  JJ = 1,ND
      COSS = COSS - DCST
      WGT = RESLT(JJ,IPOS)
      CALL HFILL(99+MODEFT,COSS,0,WGT)
      WGT = RESLT(JJ,MODEFT+4)
      CALL HFILL(109+MODEFT,COSS,0,WGT)
2132  CONTINUE
C
  250 WRITE(6,602) NPASS
      WRITE(6,605) A,SIGMAA
  605 FORMAT('0  **** A ***',1P8E12.4,/'       +/-   ',8E12.4)
C
        DO 260 K=1,8
  260   BC(K,MODEFT-1) = A(K)
  300 CONTINUE
C
        DO 220 N=1,ND
        WRITE(6,606) (RESLT(N,K),K=1,7)
  606   FORMAT(' RESLT:',F8.2,6F10.3)
  220   CONTINUE
C
C ************ ALL FITS ARE FINISHED *********
C
C    PRINT THE FINAL RESULT
C
C      WRITE(6,2271)
2271   FORMAT(1H1)
C
        DO 1000 N = 1,ND
        CHISQ = 0.
C-----   SET PARAMETERS
          DO 1050 IP=1,2
          MODEFT = IP+1
          A(IP) = FUNCTN(RESLT(1,1),N,BC(1,IP))
 1050     CONTINUE
C       A(3) = 0.55
        MODEFT = 1
        COST = RESLT(N,1)
        IXE = 14
        IF(IENDCP.EQ.1) IXE = 12
          DO 1100 IX=1,IXE
          YFIT(IX) = FUNCTN(X,IX,A)
          CHISQ = CHISQ+((SD(IX,N)-YFIT(IX))/SER(IX,N))**2
 1100     CONTINUE
C
C GRAPHICAL DISPLAY OF INPUT AND FIT
C
C     CALL HISTOG(14,SD(1,N),YFIT,0.,1.,2,4,DUMMY)
C
          DO 2491  IX = 1,IXE
          WGT = YFIT(IX)
          CALL HFILL(60+N,X(IX),0,WGT)
2491  CONTINUE
          CALL HFILL(121,COST,0,CHISQ)
C
        WRITE(6,607) COST,CHISQ,A(1),A(2),
     1           (SD(K,N),K=1,IXE),YFIT
  607   FORMAT('0COST=',F6.3,3X,'CHISQ=',E10.3,' A1-A2=',
     1        2F10.4,/' SD ',10F7.2,/,' ',4F7.2,
     2        /' YFIT ',10F7.2,/,' ',4F7.2)
 1000   CONTINUE
      CALL HISTDO
      STOP
      END
C   26/06/78 704180112  MEMBER NAME  @CERTST4 (S)           FORTRAN
C     MAIN TO TEST CERENKOV LIGHT TRACE SUBROUTINE
C---- S.YAMADA   16-11-77
C
C  GEP VERSION   4.4.87
C
C
C    @CERTST1 ON F22YAM.WORKS
C      HERE MODIFED TO REASONABLE PRINTOUT
C
C---- A MUON PENETRATES THE LEAD GLASS BLOCK.
C
      COMMON /CONSTB/ CONVMX,BK7LNG, TOTL
      COMMON /CCERSR/ SOURCE(14,25),CRNKSC(14,25),ISPCT(14,40),
     1             NEP,NACP
      DIMENSION KARRAY(1262), DIR(3)
       EQUIVALENCE (KARRAY(1),SOURCE(1,1))
C
C
      XLAMLW = .3
      XLAMHG = .7
      DXLAM = (XLAMHG-XLAMLW)/40.
C
      CALL ACORE(20000)
      CALL HBOOK1( 1,'NR VS LAMBDA,DIR 1 DEPTH 1 X0$',40,.3,.7)
      CALL HBOOK1( 2,'NR VS LAMBDA,DIR 1 DEPTH 2 X0$',40,.3,.7)
      CALL HBOOK1( 3,'NR VS LAMBDA,DIR 1 DEPTH 3 X0$',40,.3,.7)
      CALL HBOOK1( 4,'NR VS LAMBDA,DIR 1 DEPTH 4 X0$',40,.3,.7)
      CALL HBOOK1( 5,'NR VS LAMBDA,DIR 1 DEPTH 5 X0$',40,.3,.7)
      CALL HBOOK1( 6,'NR VS LAMBDA,DIR 1 DEPTH 6 X0$',40,.3,.7)
      CALL HBOOK1( 7,'NR VS LAMBDA,DIR 1 DEPTH 7 X0$',40,.3,.7)
      CALL HBOOK1( 8,'NR VS LAMBDA,DIR 1 DEPTH 8 X0$',40,.3,.7)
      CALL HBOOK1( 9,'NR VS LAMBDA,DIR 1 DEPTH 9 X0$',40,.3,.7)
      CALL HBOOK1(10,'NR VS LAMBDA,DIR 1 DEPTH 10 X0$',40,.3,.7)
      CALL HBOOK1(11,'NR VS LAMBDA,DIR 1 DEPTH 11 X0$',40,.3,.7)
      CALL HBOOK1(12,'NR VS LAMBDA,DIR 1 DEPTH 12 X0$',40,.3,.7)
      CALL HBOOK1(13,'NR VS LAMBDA,DIR 1 DEPTH 13 X0$',40,.3,.7)
      CALL HBOOK1(14,'NR VS LAMBDA,DIR 1 DEPTH 14 X0$',40,.3,.7)
C
      CALL HBOOK1(21,'NR VS LAMBDA,DIR 5 DEPTH 1 X0$',40,.3,.7)
      CALL HBOOK1(22,'NR VS LAMBDA,DIR 5 DEPTH 2 X0$',40,.3,.7)
      CALL HBOOK1(23,'NR VS LAMBDA,DIR 5 DEPTH 3 X0$',40,.3,.7)
      CALL HBOOK1(24,'NR VS LAMBDA,DIR 5 DEPTH 4 X0$',40,.3,.7)
      CALL HBOOK1(25,'NR VS LAMBDA,DIR 5 DEPTH 5 X0$',40,.3,.7)
      CALL HBOOK1(26,'NR VS LAMBDA,DIR 5 DEPTH 6 X0$',40,.3,.7)
      CALL HBOOK1(27,'NR VS LAMBDA,DIR 5 DEPTH 7 X0$',40,.3,.7)
      CALL HBOOK1(28,'NR VS LAMBDA,DIR 5 DEPTH 8 X0$',40,.3,.7)
      CALL HBOOK1(29,'NR VS LAMBDA,DIR 5 DEPTH 9 X0$',40,.3,.7)
      CALL HBOOK1(30,'NR VS LAMBDA,DIR 5 DEPTH 10 X0$',40,.3,.7)
      CALL HBOOK1(31,'NR VS LAMBDA,DIR 5 DEPTH 11 X0$',40,.3,.7)
      CALL HBOOK1(32,'NR VS LAMBDA,DIR 5 DEPTH 12 X0$',40,.3,.7)
      CALL HBOOK1(33,'NR VS LAMBDA,DIR 5 DEPTH 13 X0$',40,.3,.7)
      CALL HBOOK1(34,'NR VS LAMBDA,DIR 5 DEPTH 14 X0$',40,.3,.7)
C
      CALL HBOOK1(41,'NR VS LAMBDA,DIR 10 DEPTH 1 X0$',40,.3,.7)
      CALL HBOOK1(42,'NR VS LAMBDA,DIR 10 DEPTH 2 X0$',40,.3,.7)
      CALL HBOOK1(43,'NR VS LAMBDA,DIR 10 DEPTH 3 X0$',40,.3,.7)
      CALL HBOOK1(44,'NR VS LAMBDA,DIR 10 DEPTH 4 X0$',40,.3,.7)
      CALL HBOOK1(45,'NR VS LAMBDA,DIR 10 DEPTH 5 X0$',40,.3,.7)
      CALL HBOOK1(46,'NR VS LAMBDA,DIR 10 DEPTH 6 X0$',40,.3,.7)
      CALL HBOOK1(47,'NR VS LAMBDA,DIR 10 DEPTH 7 X0$',40,.3,.7)
      CALL HBOOK1(48,'NR VS LAMBDA,DIR 10 DEPTH 8 X0$',40,.3,.7)
      CALL HBOOK1(49,'NR VS LAMBDA,DIR 10 DEPTH 9 X0$',40,.3,.7)
      CALL HBOOK1(50,'NR VS LAMBDA,DIR 10 DEPTH 10 X0$',40,.3,.7)
      CALL HBOOK1(51,'NR VS LAMBDA,DIR 10 DEPTH 11 X0$',40,.3,.7)
      CALL HBOOK1(52,'NR VS LAMBDA,DIR 10 DEPTH 12 X0$',40,.3,.7)
      CALL HBOOK1(53,'NR VS LAMBDA,DIR 10 DEPTH 13 X0$',40,.3,.7)
      CALL HBOOK1(54,'NR VS LAMBDA,DIR 10 DEPTH 14 X0$',40,.3,.7)
C
C---- SET CONSTANTS.
      IENDCP = 0
C---- MAXIMUM CONVERSION EFFICIENCY, GOOD FOR THE NEW CONVEF
      CONVMX = 0.28
C---- TOTAL LENGTH OF THE GLASS
      TOTL = 30.0
      IF(IENDCP.EQ.1) TOTL = 25.0
      TOTLRD = TOTL/2.23
C---- BK7LNG=LIGHT GUIDE LENGH IN CM.
      BK7LNG = 6.0
      IF(IENDCP.EQ.1) BK7LNG = 50.0
      WRITE(6,611) IENDCP,TOTL,BK7LNG
  611 FORMAT('0 IENDCP,TOTL,BK7LNG=',I6,2F10.2,/)
C
C---- SET MUON  PARAMETERS
      EE = 1000.0
      T = 0.05
      JEND = TOTLRD/T
C---- CHANGE TRACK DIRECTION
      DIR(1) = 0.
      DO 1000 IDIR=1,21
      DIR(3) = 1.1-0.1*IDIR
      IF(DIR(3).GT.1.) DIR(3) = 1.
      IF(DIR(3).LT.-1.) DIR(3) = -1.
      DIR(2) = SQRT(1.-DIR(3)**2)
      WRITE(6,610) DIR
  610 FORMAT('1********* INPUT DIRECTION COSIGNS *****   ',3F8.4)
C
      SHIFT = 0.5*T*(1.-DIR(3))
C
C---- CLEAR ARRAY
        DO 400 K=1,1262
  400   KARRAY(K) = 0
C
      NTRY = 20
      TRY = NTRY
C
        DO 1 N=1,NTRY
          DO 2 J=1,JEND
          DEPTH = T*(J-1)+SHIFT
          CALL CERNKV(DEPTH,DIR,EE,T)
    2     CONTINUE
    1   CONTINUE
C
C---- HISTOGRAM OUTPUT
C     WRITE(6,600) ((SOURCE(I,J),I=1,14),J=1,25)
      RNACP = NACP
      RNACP = RNACP/TRY
      WRITE(6,603) NEP,NACP,RNACP
  603 FORMAT(' TOTAL EMITTED PHOTON  NO.=',I8,
     1     /' TOTAL ACCEPTED PHO TOELECTRON=',I8,
     2     /' PHOTO-ELE./TRIAL=',F10.2)
      WRITE(6,600) (CRNKSC(I,18),I=1,14)
  600 FORMAT('0',(10F7.0))
C---- GET PHOTO-ELECTRON NUMBERS/1 RAD.LENGTH SOURCE TRACK.
        DO 10 I=1,14
          DO 20 J=1,25
          S = SOURCE(I,J)
          IF(S.EQ.0.) S = 1.
          CRNKSC(I,J) = CRNKSC(I,J)*CONVMX/S
   20     CONTINUE
   10   CONTINUE
      WRITE(6,604)
  604 FORMAT(' PHOTO-ELECTRON NUM.PER ONE RAD.LENGTH SOURCE TRACK')
      WRITE(6,605) (CRNKSC(I,18),I=1,14)
  605 FORMAT('0',(10F7.2))
      WRITE(6,601)
  601 FORMAT(' PHOTON WAVE LENGTH DISTRIBUTION')
      WRITE(6,602) ((ISPCT(I,J),I=1,14),J=1,40)
  602 FORMAT(' ',(10I7,2X,4I7))
C
      IF(IDIR.NE.1.AND.IDIR.NE.5.AND.IDIR.NE.10) GO TO 950
      IBO = 0
      IF(IDIR.EQ.5) IBO = 20
      IF(IDIR.EQ.10) IBO = 40
      DO 940 J = 1,14
      IBO = IBO + 1
      XLAM = XLAMLW - .5*DXLAM
      DO 941  I = 1,40
      XLAM = XLAM + DXLAM
      WEIG = FLOAT(ISPCT(J,I))
      CALL HFILL(IBO,XLAM,0,WEIG)
941   CONTINUE
940   CONTINUE
C
950   CONTINUE
C
 1000 CONTINUE
      CALL HISTDO
      STOP
      END
C   26/06/78 704051735  MEMBER NAME  @CERTST6 (S)           FORTRAN
C     MAIN TO TEST CERENKOV LIGHT TRACE SUBROUTINE
C---- S.YAMADA   16-11-77
C
C   LATEST CHANGE   23.1.1986   J.OLSSON
C
C     GEP VERSION 4.4.87
C
C---- A MUON PENETRATES THE LEAD GLASS BLOCK.
C
      COMMON /CONSTB/ CONVMX,BK7LNG, TOTL
      COMMON /CCERSR/ SOURCE(20,25),CRNKSC(20,25),ISPCT(20,40),
     1             NEP,NACP
      DIMENSION KARRAY(1802), DIR(3)
       EQUIVALENCE (KARRAY(1),SOURCE(1,1))
C
      XLAMLW = .3
      XLAMHG = .7
      DXLAM = (XLAMHG-XLAMLW)/40.
C
      CALL ACORE(20000)
      CALL HBOOK1( 1,'NR VS LAMBDA,DIR 1 DEPTH 1 X0$',40,.3,.7)
      CALL HBOOK1( 2,'NR VS LAMBDA,DIR 1 DEPTH 2 X0$',40,.3,.7)
      CALL HBOOK1( 3,'NR VS LAMBDA,DIR 1 DEPTH 3 X0$',40,.3,.7)
      CALL HBOOK1( 4,'NR VS LAMBDA,DIR 1 DEPTH 4 X0$',40,.3,.7)
      CALL HBOOK1( 5,'NR VS LAMBDA,DIR 1 DEPTH 5 X0$',40,.3,.7)
      CALL HBOOK1( 6,'NR VS LAMBDA,DIR 1 DEPTH 6 X0$',40,.3,.7)
      CALL HBOOK1( 7,'NR VS LAMBDA,DIR 1 DEPTH 7 X0$',40,.3,.7)
      CALL HBOOK1( 8,'NR VS LAMBDA,DIR 1 DEPTH 8 X0$',40,.3,.7)
      CALL HBOOK1( 9,'NR VS LAMBDA,DIR 1 DEPTH 9 X0$',40,.3,.7)
      CALL HBOOK1(10,'NR VS LAMBDA,DIR 1 DEPTH 10 X0$',40,.3,.7)
      CALL HBOOK1(11,'NR VS LAMBDA,DIR 1 DEPTH 11 X0$',40,.3,.7)
      CALL HBOOK1(12,'NR VS LAMBDA,DIR 1 DEPTH 12 X0$',40,.3,.7)
      CALL HBOOK1(13,'NR VS LAMBDA,DIR 1 DEPTH 13 X0$',40,.3,.7)
      CALL HBOOK1(14,'NR VS LAMBDA,DIR 1 DEPTH 14 X0$',40,.3,.7)
      CALL HBOOK1(15,'NR VS LAMBDA,DIR 1 DEPTH 15 X0$',40,.3,.7)
      CALL HBOOK1(16,'NR VS LAMBDA,DIR 1 DEPTH 16 X0$',40,.3,.7)
      CALL HBOOK1(17,'NR VS LAMBDA,DIR 1 DEPTH 17 X0$',40,.3,.7)
      CALL HBOOK1(18,'NR VS LAMBDA,DIR 1 DEPTH 18 X0$',40,.3,.7)
      CALL HBOOK1(19,'NR VS LAMBDA,DIR 1 DEPTH 19 X0$',40,.3,.7)
      CALL HBOOK1(20,'NR VS LAMBDA,DIR 1 DEPTH 20 X0$',40,.3,.7)
C
      CALL HBOOK1(21,'NR VS LAMBDA,DIR 5 DEPTH 1 X0$',40,.3,.7)
      CALL HBOOK1(22,'NR VS LAMBDA,DIR 5 DEPTH 2 X0$',40,.3,.7)
      CALL HBOOK1(23,'NR VS LAMBDA,DIR 5 DEPTH 3 X0$',40,.3,.7)
      CALL HBOOK1(24,'NR VS LAMBDA,DIR 5 DEPTH 4 X0$',40,.3,.7)
      CALL HBOOK1(25,'NR VS LAMBDA,DIR 5 DEPTH 5 X0$',40,.3,.7)
      CALL HBOOK1(26,'NR VS LAMBDA,DIR 5 DEPTH 6 X0$',40,.3,.7)
      CALL HBOOK1(27,'NR VS LAMBDA,DIR 5 DEPTH 7 X0$',40,.3,.7)
      CALL HBOOK1(28,'NR VS LAMBDA,DIR 5 DEPTH 8 X0$',40,.3,.7)
      CALL HBOOK1(29,'NR VS LAMBDA,DIR 5 DEPTH 9 X0$',40,.3,.7)
      CALL HBOOK1(30,'NR VS LAMBDA,DIR 5 DEPTH 10 X0$',40,.3,.7)
      CALL HBOOK1(31,'NR VS LAMBDA,DIR 5 DEPTH 11 X0$',40,.3,.7)
      CALL HBOOK1(32,'NR VS LAMBDA,DIR 5 DEPTH 12 X0$',40,.3,.7)
      CALL HBOOK1(33,'NR VS LAMBDA,DIR 5 DEPTH 13 X0$',40,.3,.7)
      CALL HBOOK1(34,'NR VS LAMBDA,DIR 5 DEPTH 14 X0$',40,.3,.7)
      CALL HBOOK1(35,'NR VS LAMBDA,DIR 5 DEPTH 15 X0$',40,.3,.7)
      CALL HBOOK1(36,'NR VS LAMBDA,DIR 5 DEPTH 16 X0$',40,.3,.7)
      CALL HBOOK1(37,'NR VS LAMBDA,DIR 5 DEPTH 17 X0$',40,.3,.7)
      CALL HBOOK1(38,'NR VS LAMBDA,DIR 5 DEPTH 18 X0$',40,.3,.7)
      CALL HBOOK1(39,'NR VS LAMBDA,DIR 5 DEPTH 19 X0$',40,.3,.7)
      CALL HBOOK1(40,'NR VS LAMBDA,DIR 5 DEPTH 20 X0$',40,.3,.7)
C
      CALL HBOOK1(41,'NR VS LAMBDA,DIR 10 DEPTH 1 X0$',40,.3,.7)
      CALL HBOOK1(42,'NR VS LAMBDA,DIR 10 DEPTH 2 X0$',40,.3,.7)
      CALL HBOOK1(43,'NR VS LAMBDA,DIR 10 DEPTH 3 X0$',40,.3,.7)
      CALL HBOOK1(44,'NR VS LAMBDA,DIR 10 DEPTH 4 X0$',40,.3,.7)
      CALL HBOOK1(45,'NR VS LAMBDA,DIR 10 DEPTH 5 X0$',40,.3,.7)
      CALL HBOOK1(46,'NR VS LAMBDA,DIR 10 DEPTH 6 X0$',40,.3,.7)
      CALL HBOOK1(47,'NR VS LAMBDA,DIR 10 DEPTH 7 X0$',40,.3,.7)
      CALL HBOOK1(48,'NR VS LAMBDA,DIR 10 DEPTH 8 X0$',40,.3,.7)
      CALL HBOOK1(49,'NR VS LAMBDA,DIR 10 DEPTH 9 X0$',40,.3,.7)
      CALL HBOOK1(50,'NR VS LAMBDA,DIR 10 DEPTH 10 X0$',40,.3,.7)
      CALL HBOOK1(51,'NR VS LAMBDA,DIR 10 DEPTH 11 X0$',40,.3,.7)
      CALL HBOOK1(52,'NR VS LAMBDA,DIR 10 DEPTH 12 X0$',40,.3,.7)
      CALL HBOOK1(53,'NR VS LAMBDA,DIR 10 DEPTH 13 X0$',40,.3,.7)
      CALL HBOOK1(54,'NR VS LAMBDA,DIR 10 DEPTH 14 X0$',40,.3,.7)
      CALL HBOOK1(55,'NR VS LAMBDA,DIR 10 DEPTH 15 X0$',40,.3,.7)
      CALL HBOOK1(56,'NR VS LAMBDA,DIR 10 DEPTH 16 X0$',40,.3,.7)
      CALL HBOOK1(57,'NR VS LAMBDA,DIR 10 DEPTH 17 X0$',40,.3,.7)
      CALL HBOOK1(58,'NR VS LAMBDA,DIR 10 DEPTH 18 X0$',40,.3,.7)
      CALL HBOOK1(59,'NR VS LAMBDA,DIR 10 DEPTH 19 X0$',40,.3,.7)
      CALL HBOOK1(60,'NR VS LAMBDA,DIR 10 DEPTH 20 X0$',40,.3,.7)
C
C---- SET CONSTANTS.
      IENDCP = 0
C---- MAXIMUM CONVERSION EFFICIENCY, GOOD FOR THE NEW CONVEF
      CONVMX = 0.28
C---- TOTAL LENGTH OF THE GLASS
      TOTL = 30.0
      IF(IENDCP.EQ.1) TOTL = 25.0
C     TOTLRD = TOTL/2.23
      TOTLRD = TOTL/1.53
C---- BK7LNG=LIGHT GUIDE LENGH IN CM.
C     BK7LNG = 6.0
      BK7LNG = 0.0
      IF(IENDCP.EQ.1) BK7LNG = 50.0
      WRITE(6,611) IENDCP,TOTL,BK7LNG
  611 FORMAT('  IENDCP,TOTL,BK7LNG=',I6,2F10.2,/)
C
C---- SET MUON  PARAMETERS
      EE = 1000.0
      T = 0.05
      JEND = TOTLRD/T
C---- CHANGE TRACK DIRECTION
      DIR(1) = 0.
C     DO 1000 IDIR=1,5
      DO 1000 IDIR=1,21
C     DIR(3) = 1.2-0.2*IDIR
C     DIR(3) = 1.05-0.05*IDIR
      DIR(3) = 1.10-0.10*IDIR
      IF(DIR(3).GT.1.) DIR(3) = 1.
      IF(DIR(3).LT.-1.) DIR(3) = -1.
      DIR(2) = SQRT(1.-DIR(3)**2)
      WRITE(6,610) DIR
  610 FORMAT('1********* INPUT DIRECTION COSIGNS *****   ',3F8.4)
C
      SHIFT = 0.5*T*(1.-DIR(3))
C
C---- CLEAR ARRAY
        DO 400 K=1,1802
  400   KARRAY(K) = 0
C
      NTRY = 100
C     NTRY = 20
      TRY = NTRY
C
        DO 1 N=1,NTRY
          DO 2 J=1,JEND
          DEPTH = T*(J-1)+SHIFT
          CALL CERNKV(DEPTH,DIR,EE,T)
    2     CONTINUE
    1   CONTINUE
C
C---- HISTOGRAM OUTPUT
C     WRITE(6,600) (SOURCE(I,18),I=1,20)
      RNACP = NACP
      RNACP = RNACP/TRY
      WRITE(6,603) NEP,NACP,RNACP
  603 FORMAT(' TOTAL EMITTED PHOTON  NO.=',I8,
     1     /' TOTAL ACCEPTED PHO TOELECTRON=',I8,
     2     /' PHOTO-ELE./TRIAL=',F10.2)
C//////   RESULTS TO BE USED IN FITTING PROGRAM   ///////////////
      WRITE(6,1603) DIR(3)
 1603 FORMAT(' ',F6.3)
      WRITE(6,600) (CRNKSC(I,18),I=1,20)
  600 FORMAT(' ',(10F7.0))
C///////////////////////////////////////////////////////////////
C---- GET PHOTO-ELECTRON NUMBERS/1 RAD.LENGTH SOURCE TRACK.
        DO 10 I=1,20
          DO 20 J=1,25
          S = SOURCE(I,J)
          IF(S.EQ.0.) S = 1.
          CRNKSC(I,J) = CRNKSC(I,J)*CONVMX/S
   20     CONTINUE
   10   CONTINUE
      WRITE(6,604)
  604 FORMAT(' PHOTO-ELECTRON NUM.PER ONE RAD.LENGTH SOURCE TRACK')
      WRITE(6,605) (CRNKSC(I,18),I=1,20)
  605 FORMAT(' ',(10F7.2))
      WRITE(6,601)
  601 FORMAT(' PHOTON WAVE LENGTH DISTRIBUTION')
      WRITE(6,602) ((ISPCT(I,J),I=1,20),J=1,40)
C 602 FORMAT(' ',(10I7,2X,4I7))
  602 FORMAT(' ',10I7)
C
      IF(IDIR.NE.1.AND.IDIR.NE.5.AND.IDIR.NE.10) GO TO 950
      IBO = 0
      IF(IDIR.EQ.5) IBO = 20
      IF(IDIR.EQ.10) IBO = 40
      DO 940 J = 1,20
      IBO = IBO + 1
      XLAM = XLAMLW - .5*DXLAM
      DO 941  I = 1,40
      XLAM = XLAM + DXLAM
      WEIG = FLOAT(ISPCT(J,I))
      CALL HFILL(IBO,XLAM,0,WEIG)
941   CONTINUE
940   CONTINUE
C
950   CONTINUE
C
 1000 CONTINUE
      CALL HISTDO
      STOP
      END
C   26/06/78 601232211  MEMBER NAME  @FNPECER (S)           FORTRAN
C
C    PLOT THE FUNCTION NPECER
C         GEP VERSION
C
C
      DIMENSION X(14),A(2)
      DATA X/0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5,
     1  13.4/
C
      LOGICAL*1 TEXT
      DIMENSION TEXT(72)
C
C
      IPRO = 0
      NPTS = 14
C
C PLOT INITIALIZATION
C
      CALL ACORE(10000)
      DO 3705  I = 1,21
      CALL HBOOK1( I,' NR OF PHOTONS VS. RAD. LENGTH$',20,0.,20.)
3705  CONTINUE
      DO 3706  I = 1,2
      CALL HBOOK1(110+ I,' FIT PARAMETER, GLOBAL FIT$',21,-1.05,1.05)
3706  CONTINUE
C
7608  FORMAT(72A1)
7508  FORMAT(' ',72A1)
C     READ 7608, TEXT
C     WRITE(6,7508) TEXT
C
C  LOOP OVER VARIOUS COST STARTING HERE   ###########################
C
      COST = 1.10
      DELCO = 0.10
      ND = 0
C
 1    CONTINUE
      COST = COST - DELCO
      IF(COST.LT.-1.) GO TO 1000
      ND = ND + 1
C
      DO 1232  JJ = 1,NPTS
      CALL NPECER(1,X(JJ),COST,PENUM,A(1))
      WGT = PENUM
      CALL HFILL(ND,X(JJ),0,WGT)
1232  CONTINUE
C
      WGT = A(1)
      CALL HFILL(111,COST,0,WGT)
      WGT = A(2)
      CALL HFILL(112,COST,0,WGT)
C
      GO TO 1
C
 1000   CONTINUE
      CALL HISTDO
      STOP
      END
C   26/06/78 705112054  MEMBER NAME  @NPECER  (S)           FORTRAN
C
C    PLOT THE FUNCTION NPECER
C         GEP VERSION
C
C
      DIMENSION X(14),A(2)
      DATA X/0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5,
     1  13.4/
C
      LOGICAL*1 TEXT
      DIMENSION TEXT(72)
C
C
      IPRO = 0
      NPTS = 14
C
C PLOT INITIALIZATION
C
      CALL ACORE(10000)
      DO 3705  I = 1,21
      CALL HBOOK1( I,' NR OF PHOTONS VS. RAD. LENGTH$',20,0.,20.)
3705  CONTINUE
C
7608  FORMAT(72A1)
7508  FORMAT(' ',72A1)
C     READ 7608, TEXT
C     WRITE(6,7508) TEXT
C
C  LOOP OVER VARIOUS COST STARTING HERE   ###########################
C
      COST = 1.10
      DELCO = 0.10
      ND = 0
C
 1    CONTINUE
      COST = COST - DELCO
      IF(COST.LT.-1.) GO TO 1000
      ND = ND + 1
C
      DO 1232  JJ = 1,NPTS
      CALL NPECER(1,X(JJ),COST,PENUM)
      WGT = PENUM
      CALL HFILL(ND,X(JJ),0,WGT)
1232  CONTINUE
C
      GO TO 1
C
 1000   CONTINUE
      CALL HISTDO
      STOP
      END
