	subroutine ice
c	read in data sets for this t value
c
	common/icenot/tval(50),nddt
	common/itadd/itadd(7)
	common/lbias/lbias(7)
	include '/u1/gcf/e350/crystal/cb/pperp'
	include '/u1/gcf/e350/crystal/cb/foxin'
	include '/u1/gcf/e350/crystal/cb/datbin'
	include '/u1/gcf/e350/crystal/cb/norm'
	include '/u1/rgk/icicle/cb/ztpltc'
	include '/u1/gcf/e350/crystal/cb/clouds'
	include '/u1/gcf/e350/crystal/cb/abssiz'
	include '/u1/gcf/e350/crystal/cb/type'
	include '/u1/gcf/e350/crystal/cb/irspec'
	include '/u1/gcf/e350/crystal/cb/it'
c       equivalence(it,itr)
	include '/u1/gcf/e350/crystal/cb/parms'
	include '/u1/gcf/e350/crystal/cb/mceffy'
	dimension atemp(180),btemp(180),ctemp(180),fudfl(180),fudmt(180),
& 	dtemp(180),etemp(180),ftemp(180),ebegg(6)
	double precision avtd,fdd,vald,ppind,eintd,abd
c
c	modified b130, b170, e150 from 130,170,135 to 140,180,140 on dec11 85
        data notrgo/0/,nony/0/,ebegg/100.,140.,180.,90.,120.,140./
c
7	format(' range of fit'2f7.3' average'f8.4' events 'f9.2)
8	format(' efficiencies',12f8.4/(13x,12f8.4))
9	format('sets empty'i3' full'i3' t range'2i3' full/empty fudges'2g11.3)
100	format('energy range of data is 'f5.1' - 'f5.1 ' GeV. No data bins='i4)
101	format('energy range of fit is 'f5.1' - 'f5.1' GeV  bins'i4' -'i4' no='
& 	i3)
102	format('pperp/e   '15f6.4/(11x,15f6.4))
104	format('for trigger no'i3' at energy of 'f6.2' GeV and pperp/e of 'f6.3
& 	' the tmin correction is'f6.4/'full and empty acceptances are'2f8.4)
105	format('for trigger no 'i3' lower range of fit increased to'f10.3)
106	format('empty/full raw cross section ratio = 'f8.4)
107	format('full cross section ='f10.4' empty ='f10.4)
108	format('cannot add sets no'i2' sets are'2i2)
109	format('Stopping because it1 it2 no of t bins ='3i3)
110	format('average value of t or pperp2 ='f9.4)
140	format('Stopping on reading CACTUS file.  it1 it2 ='2i5)
141	format('**')
c
	if (nadd.gt.0.and.notrgo.ne.0) ntrdis=notrgo
	if (nadd.gt.0.and.nony.ne.0) notrig=nony
19	it1=1
	it2=-10
	fxlo=1.e-7
	fxhi=1.e10
	do 30 i=1,ntrdis
	if (e1o(i).gt.0.) e1(i)=e1o(i)
30	if (e2o(i).gt.0.) e2(i)=e2o(i)
	iwset=1
	uplabl=.false.
c
c       old namelist cactus
c	read(2,*,end=50) combin,iocon,iocon1,iocon2,it1,it2,dofit,iglow,uplabl
	read(2,*,end=50) it1,it2
	if (it1.le.0) goto 50
	do 3 i=1,ntrdis
3	ebeg(i)=ebegg(lbias(i))
	if (ntr.eq.1) jrspec=0
	if(it2.lt.0) it2=it1
	nddt1=nddt-1
	if (it1+itadd(notrig).gt.nddt1.or.it2+itadd(notrig).gt.nddt) goto 51
c
c	counters to find average t or pperp2
	absum=0.
	avtsum=0.
c
c	set up tlo and thi - note that tval is set up for the final trigger
	tlo=-tval(it1+itadd(notrig))
	thi=-tval(it2+itadd(notrig)+1)
c	set up fixed(1)
	if (nfs.or.-thi.le.0.41) then
	fixed(1)=.true.
	else
	fixed(1)=.false.
	endif
c
c	set up usei2
	usei2=.false.
c	if (.not.nfs.and.ipart.eq.1.and.-thi.lt.0.51) usei2=.true.
c
c	set up ntr
	ntr=1
c	if (.not.nfs.and.-thi.lt.0.31) ntr=2
c
	if (nfs.and.-tlo.gt.0.24) combin(1)=2
	print 300,it1,it2,combin,iocon,iocon1(1),iocon2(1),iglow,(ebeg(i),
&	i=1,ntrdis)
300	format('CACTUS it1 it2'2i3' combin'7i2' iocon'7i2' iocon1 iocon2'
&	2i2' iglow'l2' ebeg'7f6.1)
c
	do 1 itr=1,ntrdis
c	read in target full data
	call loadit(dtemp,etemp,ftemp,elow(itr),de(itr),.false.,
&	combin(itr),core(itr),noe(itr),fudfl,itadd(itr))
	e2(itr)=amin1(e2(itr),elow(itr)+float(noe(itr))*de(itr))
	ao=0.
	ppin=pinc*core(itr)
	nume=noe(itr)
	if (.not.mt(itr)) goto 2
	iwseto=iwset
	iwset=iwset+1
	print 9,iwset,iwseto,it1,it2,fulinv(itr),emtinv(itr)
c	read in empty target data
	call loadit(atemp,btemp,ctemp,dum1,dum2,.true.,combin(itr),
&	corem(itr),nume,fudmt,itadd(itr))
	fuladd=0.
	empadd=0.
2	idoon=0
	do 6 i=1,nume
	eint=elow(itr)+de(itr)*(i-0.5)
	if (.not.mt(itr)) goto 17
c	see if mt or full acceptance too low
	if (fudfl(i).lt.0..or.fudmt(i).lt.0.) goto 12
	a1=dtemp(i)*fulinv(itr)
	a2=-atemp(i)*emtinv(itr)
	fuladd=fuladd+a1
	empadd=empadd-a2
	a=a1+a2
	if (a1.ne.a2) ftemp(i)=(a1*ftemp(i)-a2*ctemp(i))/(a1-a2)
	dtemp(i)=a
	et=sqrt(abs(ful2(itr)*etemp(i)+emt2(itr)*btemp(i)))
	fdiv=correc(eint*core(itr),iacc(itr),tlo,thi)
	if (fdiv.lt.1.e-10) fdiv=fudfl(i)
	ftest(i,itr)=etest(itr)/fdiv
	if (et.gt.0.) goto 36
c	need to calculate acceptance when error is zero
	etemp(i)=1./ftest(i,itr)
	goto 35
36	etemp(i)=1./et
35	ao=ao+dtemp(i)
	ptoe(i,itr)=ftemp(i)
	binb(i,itr)=dtemp(i)
	ebinb(i,itr)=etemp(i)
	if (ftemp(i).lt.-99..and.eint.gt.e1(itr).and.eint.lt.e2(itr)) goto
& 	31
	goto 6
c	for case when no empty subtraction
17	if (fudfl(i).lt.0.) goto 12
	dtemp(i)=dtemp(i)*fulinv(itr)
	et=sqrt(ful2(itr)*etemp(i))
	ftest(i,itr)=etest(itr)/correc(eint*core(itr),iacc(itr),tlo,thi)
	if (et.gt.0.) goto 36
	etemp(i)=1./ftest(i,itr)
	goto 35
c
c	acceptance is too low
12	print 104,itr,eint,ftemp(i),fudstr,fudfl(i),fudmt(i)
	if (eint.lt.e1(itr).or.(eint.gt.e2(itr).and.eint.ge.90.)) goto 24
31	if (eint-e1(itr).le.10.0) goto 37
16	if (.not.cross) goto 80
37	e1(itr)=eint+0.5*de(itr)
	ebeg(itr)=amax1(e1(itr),ebeg(itr))
	idoon=1
	if (nadd.eq.0) goto 24
	if (itr.eq.iadb(1)) e1(iada(1))=e1(iadb(1))
	if (nadd.ge.2.and.itr.eq.iadb(2)) e1(iada(2))=e1(iadb(2))
	if (itr.eq.iadb(1)) ebeg(iada(1))=ebeg(iadb(1))
	if (nadd.ge.2.and.itr.eq.iadb(2)) ebeg(iada(2))=ebeg(iadb(2))
24	binb(i,itr)=-1.e-5
	ebinb(i,itr)=0.
6	continue
	if (idoon.eq.1) print 105,itr,e1(itr)
	if (.not.mt(itr)) goto 10
	if (fuladd.le.0.) goto 18
	divide=empadd/fuladd
	print 106,divide
	goto 10
18	print 107,fuladd,empadd
c
c	check efficiencies are large enuf to permit fitting of points
10	if (.not.dofit) goto 23
	do 20 i=1,nume
	if (elow(itr)+float(i)*de(itr).lt.e1(itr)) goto 20
	if (abs(binb(i,itr)+1.).gt.1.e-6) goto 21
	e1(itr)=elow(itr)+float(i)*de(itr)
20	continue
21	do 22 i=1,nume
	j=nume-i+1
	if (elow(itr)+float(j)*de(itr).gt.e2(itr)) goto 22
	if (abs(binb(j,itr)+1.).gt.1.e-6) goto 23
	e2(itr)=elow(itr)+float(j-1)*de(itr)
22	continue
23	if (e1(itr).lt.elow(itr)) e1(itr)=elow(itr)
	ie1(itr)=(e1(itr)-elow(itr)+0.99*de(itr))/de(itr)
	ie1(itr)=ie1(itr)+1
	e1(itr)=elow(itr)+de(itr)*(ie1(itr)-1)
	ehgh=elow(itr)+noe(itr)*de(itr)
	if (e2(itr).gt.ehgh) e2(itr)=ehgh
	ie2(itr)=(e2(itr)-elow(itr))/de(itr)
	ie2(itr)=min0(nume,ie2(itr))
	e2(itr)=elow(itr)+de(itr)*ie2(itr)
	moe(itr)=ie2(itr)-ie1(itr)+1
	if (moe(itr).le.0.and.dofit.and.(.not.cross).and.itr.le.notrig)
& 	call abort
	ie1(itr)=min0(ie1(itr),ie2(itr))
	print 100,elow(itr),ehgh,noe(itr)
	print 101,e1(itr),e2(itr),ie1(itr),ie2(itr),moe(itr)
c	calculate average t values for bins being fitted
	ia=ie1(itr)
	ib=ie2(itr)
	avtd=0.
	abd=0.
	if (moe(itr).le.0) goto 28
	do 15 ie=ia,ib
	vald=binb(ie,itr)
	if (vald.le.1.e-5) goto 15
	abd=abd+vald
c	ftemp(ie) is same as ptoe(ie,itr)
	eintd=(elow(itr)+(float(ie)-.5)*de(itr))*core(itr)
	if (dopp) avtd=avtd+(ftemp(ie)*eintd)**2*vald
	ppind=pinc*core(itr)
	fdd=ftemp(ie)**2
	if (.not.dopp) avtd=avtd+(am2id+am2jd-2.*dsqrt(ppind**2+am2jd)*eintd
& 	+2.*ppind*dsqrt(eintd**2*(1.-fdd)-am2id))*vald
15	continue
	avtsum=avtsum+avtd
	absum=absum+abd
	if (abd.gt.0.) avt=avtd/abd
28	if (abd.le.0.) avt=0.5*(tlo+thi)
	print 7,tlo,thi,avt,ao
	print 102,(ptoe(i,itr),i=1,nume)
	iwset=iwset+1
1	continue
        it=itr
	if (absum.gt.0.) avt=avtsum/absum
	if(absum.le.0 .and. dopp) call abort
	if (absum.le.0.) avt=.5*(tlo+thi)
	print 110,avt
	if (nadd.eq.0) return
	if (notrgo.ne.0) goto 25
	do 32 k=1,nadd
	if (iada(k).gt.ntrdis) goto 33
	if (iadb(k).gt.ntrdis) goto 33
32	continue
	goto 34
33	print 108,k,iada(k),iadb(k)
	nadd=0
	return
34	notrgo=ntrdis
	nuts=notrgo-nadd
	if (lwbias(1).gt.iadb(1)) lwbias(1)=lwbias(1)-1
	if (lwbias(2).gt.iadb(1)) lwbias(2)=lwbias(2)-1
	if (nadd.ge.2.and.lwbias(1).gt.iadb(2)) lwbias(1)=lwbias(1)-1
	if (nadd.ge.2.and.lwbias(2).gt.iadb(2)) lwbias(2)=lwbias(2)-1
	do 29 k=1,nadd
	icb=iadb(k)
	icj=notrgo-k
	do 29 j=icb,icj
29	igrp(j)=igrp(j+1)
25	ntrdis=nuts
	nony=notrig
	do 40 k=1,nadd
	if (iada(k).le.nony) notrig=notrig-1
40	continue
	do 27 k=1,nadd
	ica=iada(k)
	icb=iadb(k)
	num=noe(ica)
	icj=notrgo-k
	do 26 j=icb,icj
	ie1(j)=ie1(j+1)
	ie2(j)=ie2(j+1)
	elow(j)=elow(j+1)
	de(j)=de(j+1)
	noe(j)=noe(j+1)
	e1(j)=e1(j+1)
	e2(j)=e2(j+1)
	ebeg(j)=ebeg(j+1)
26	moe(j)=moe(j+1)
	do 11 i=1,num
	if (ebinb(i,ica).eq.0.or.ebinb(i,icb).eq.0) goto 14
	ebinb(i,ica)=1./sqrt(1./ebinb(i,ica)**2+1./ebinb(i,icb)**2)
	if (binb(i,ica).le.0.) goto 4
	if (binb(i,icb).le.0.) goto 11
	ptoe(i,ica)=(binb(i,ica)*ptoe(i,ica)+binb(i,icb)*ptoe(i,icb))/
& 	(binb(i,ica)+binb(i,icb))
5	binb(i,ica)=binb(i,ica)+binb(i,icb)
	goto 11
4	if (binb(i,icb).le.0.) goto 5
	binb(i,ica)=binb(i,icb)
	ptoe(i,ica)=ptoe(i,icb)
	goto 11
14	ebinb(i,ica)=0.
	binb(i,ica)=-1.e-5
	ptoe(i,ica)=-100.
11	continue
	do 13 j=icb,icj
	num=noe(j)
	do 13 i=1,num
	binb(i,j)=binb(i,j+1)
	ebinb(i,j)=ebinb(i,j+1)
	ptoe(i,j)=ptoe(i,j+1)
	ftest(i,j)=ftest(i,j+1)
13	continue
27	continue
	return
c
80      call abort
c
c	it1 or it2 outside t range read in so print message and stop'
51	print 109,it1,it2,nddt1
c
c	ending
50	close(2)
cnono	if (dolist) call hecuba
	if (doplot) write(3,141)
	print 140,it1,it2
	stop
	end
