	subroutine setpar
c	sets the initial parameters
c
	include '/u1/gcf/e350/crystal/cb/irspec'
	include '/u1/gcf/e350/crystal/cb/hotsun'
	include '/u1/gcf/e350/crystal/cb/it'
	include '/u1/gcf/e350/crystal/cb/igr'
	include '/u1/gcf/e350/crystal/cb/type'
	include '/u1/gcf/e350/crystal/cb/foxin'
	include '/u1/gcf/e350/crystal/cb/datbin'
	include '/u1/gcf/e350/crystal/cb/parms'
	include '/u1/gcf/e350/crystal/cb/pperp'
	include '/u1/gcf/e350/crystal/cb/abssiz'
	include '/u1/gcf/e350/crystal/cb/pperpd'
	include '/u1/gcf/e350/crystal/cb/clouds'
	include '/u1/gcf/e350/crystal/cb/x'
	include '/u1/gcf/e350/crystal/cb/where'
	include '/u1/gcf/e350/crystal/cb/wenbin'
	include '/u1/gcf/e350/crystal/cb/beam'
	include '/u1/gcf/e350/crystal/cb/const'
	external regge,resp,confit,bmfunc,tfixed,cspi,cset
	logical fitbmo,ftsigo,ftalpo
c
100	format(' doing preliminary fit with alpha fixed at 'f8.4)
101	format(' alpha for f term fixed at 'f8.4)
102	format(' fixed mass normalization fixed at 'f7.1)
103	format(' f normalization = '2f13.4)
104	format(' alpha fixed at 'f9.4' ipart avt'i3,f9.4)
c
	s=2.*amp*sqrt(pinc**2+am2j)+amp2+am2j
	if (.not.fixed(1)) goto 5
c	use empirically determined cross-section from fitting 70<e<102
	if (.not.dopp) goto 10
c	set tlp and thp as t range of bin (negative) at x=1
	eee=pinc**2+am2j
	tlp=am2i+am2j-2.*eee+2.*pinc*sqrt(eee+tlo-am2i)
	thp=am2i+am2j-2.*eee+2.*pinc*sqrt(eee+thi-am2i)
	goto 50
10	tlp=tlo
	thp=thi
50	if (ipart.eq.1) then
		if (nfs) then
		pff=xsec(s,tlp,thp,cspi,1)/(fdmicr*bratio(1))
		fxlo=pff
		fxhi=4.*pff
	print 888,fxlo,fxhi
888	format(' fxlo fxhi'2f9.3)
		else
		pff=xsec(s,tlp,thp,cspi,2)/(fdmicr*bratio(1))
		endif
	else
		if (nfs) then
		pff=xsec(s,tlp,thp,cset,3)/(fdmicr*bratio(2))
		else
		pff=xsec(s,tlp,thp,cset,4)/(fdmicr*bratio(2))
		endif
	endif
51	print 102,pff
5	moint=1
	do 2 i=1,notrig
2	call wgtint(0.,de(i),moint,val1(1,i),wgt1(1,i))
	nointr=8
	noinfx=1
	nobeam=1
	if (fitbem) nobeam=3
	if (nfs) goto 31
	if (abs(e1(1)-140.).lt.1.) goto 18
	if (usei2.or.(.not.usei2.and.abs(thi).gt.0.51)) alfain(1)=
& 	-0.448+0.718*exp(-1.19*avt-5.8*avt**2)
	if (.not.usei2.and.abs(thi).lt.0.51) alfain(1)=-0.460+0.716*
& 	exp(0.06*avt-2.92*avt**2)
c	use 100 gev fit
	alfain(1)=-0.448+0.718*exp(-1.19*avt-5.8*avt**2)
c	200 gev fit over entire range
c	alfain(1)=-0.369+0.754*exp(-.27*avt-2.29*avt**2)+.052*avt
	goto 17
18	alfain(1)=-0.746+1.143*exp(0.784*avt-.421*avt**2)
c	use 100 GeV curve
	alfain(1)=-0.721+1.011*exp(0.62*avt-1.00*avt**2)
c
c	200 gev fit
	alfain(1)=-0.405+0.726*exp(-.19*avt-3.14*avt**2)+0.0880*avt
	goto 17
31	if (ipart.eq.2) goto 41
	if (e1(1).lt.60.) goto 32
c	use 100 gev curve
	alfain(1)=-0.691+1.29*exp(1.40*avt-0.22*avt**2)
c	fit to 200 GeV
	alfain(1)=-0.575+1.077*exp(1.73*avt-0.55*avt**2)
	goto 17
32	alfain(1)=-0.675+1.275*exp(1.277*avt-1.03*avt**2)
	goto 17
41	alfain(1)=-0.50+0.98*exp(0.49*avt-0.7*avt**2)
17	print 104,alfain(1),ipart,avt
12	continue
	alfain(2)=alfain(1)
	call bemset
	if(.not. dofit) return
	if (usei2) call subi2
c
c	if (ftbrem) goto 9
	if (jrspec.eq.1) goto 20
	call regset(e1(1),e2(1),ptoe(ie1(1),1),ptoe(ie2(1),1),1)
	x1=x1a(1)
	x2=x2a(1)
21	if (nfs.and.fitbem) goto 1
c	assume that triple regge term dominates for all inclusive reaction
c	assume first term dominant if ntr gt 1
	call sunny
	if (.not.ftnorm(1)) goto 7
	tot=0.
	dattot=0.
	alpha=alfain(1)
	const=conal(1)
	para=0.
	noint=nointr
	do 8 it=1,notrig
	ello=e1(it)
	ia=ie1(it)
	ib=ie2(it)
	do 3 i=ia,ib
	pperpd=ptoe(i,it)
	call reglim(ello,de(it),ido,0,fixed(igrp(it)))
	if (ido.eq.0) tot=tot+calct(ello,regge,resp,bmfunc)
	dattot=dattot+binb(i,it)
3	ello=ello+de(it)
	para=dattot/tot+para
	if (fixed(1).and.fitfix) param(ifixed)=pff
8	continue
	para=para/float(notrig)
	param(inorm(1))=para
	if (ntr.eq.1) return
	do 6 i=2,ntr
6	if (ftnorm(i)) param(inorm(i))=.01*para
	if (.not.nfs) goto 7
	if (jrspec.eq.0) return
	param(inorm(2))=param(inorm(1))
	alfain(2)=alfain(1)-0.2
	return
c
7	if (ftnorm(2)) param(inorm(2))=.5
	if (jpart(1).gt.2) goto 16
	if (ipart.eq.1) plim2(1)=.71*.20*1000.*enhghi((nobeam+1)/2)*xsec(
& 	s,tlo,thi,cspi,0)/(cspi(s,0.)*fdmicr)
	if (ipart.eq.2) plim2(1)=.20*1000.*enhghi((nobeam+1)/2)*xsec(
& 	s,tlo,thi,cset,0)*bratio(1)/(cspi(s,0.)*bratio(2)*fdmicr)
	if (jpart(1).eq.2) plim2(1)=0.65*plim2(1)
	if (jpart(2).eq.2) plim2(2)=0.65*plim2(1)
	print 103,plim2
	return
16	plim2(2)=0.
	plim2(1)=0.45
	if (jpart(1).lt.3) return
	amula=1.+conal(2)-2.*alfain(1)
	amulb=1.+conal(1)-2.*alfain(1)
	if (jpart(1).ge.5) goto 24
	if (ipart.eq.1) plim2(1)=0.7*exp(9.7*avt)*(x1**amula-x2**amula)*
& 	amulb/(amula*(x1**amulb-x2**amulb)*sqrt(s))
	if (ipart.eq.2) plim2(1)=0.7*exp(12.*avt)*(x1**amula-x2**amula)*
& 	amulb/(amula*(x1**amulb-x2**amulb)*sqrt(s))
25	print 103,plim2
	return
24	plim2(1)=0.9*(x1**amula-x2**amula)*amulb/(amula*(x1**amulb-
& 	x2**amulb)*sqrt(s))
	goto 25
c
1	if (.not.fitfix) goto 13
	if (ntr.gt.0) goto 13
	call sunny
	para=0.
	noint=noinfx
	do 11 it=1,notrig
11	para=para+parfix(e1(it),ie1(it),ie2(it),de(it),binb(1,it))
	param(ifixed)=para/float(notrig)
	return
c
c	first fit with alpha fixed to get the triple regge terms normalization
13	alpha=alfain(1)
	print 100,alpha
	ftalpo=fitalp(1)
	fitalp(1)=.false.
	fitbmo=fitbem
	ftsigo=fitsga
	fitbem=.false.
	fitsga=.false.
	ntrsav=ntr
	ntr=1
	call sunny
	noint=noinfx
	if (.not.fitfix) goto 30
c	para=0.
c	pfsave=pff
c	pff=1.0
c	do 4 it=1,notrig
c	igr=jpart(igrp(it))
c4	para=para+parfix(e1(it),ie1(it),ie2(it),de(it),binb(1,it))
c	param(ifixed)=para/float(notrig)
c	pff=pfsave
c	changed feb25 85  commented out above 9 lines and added following line
	param(ifixed)=pff
	goto 40
c
c	calculate triple regge normalization by subtracting fixed mass
c	term from total
30	ia=ie1(1)
	ib=ie2(1)
	totr=0.
	tot=0.
	dattot=0.
	ello=e1(1)
	it=1
	const=conal(1)
	alpha=alfain(1)
	igr=jpart(igrp(1))
	pfsave=pff
	pff=1.0
	do 35 i=ia,ib
	call fixlim(ello,de(1),ido)
	if (ido.eq.0) tot=tot+calct(ello,tfixed,resp,bmfunc)
	pperpd=ptoe(i,it)
	call reglim(ello,de(1),ido,0,fixed)
	if (ido.eq.0) totr=totr+calct(ello,regge,resp,bmfunc)
	dattot=dattot+binb(i,1)
35	ello=ello+de(1)
	param(inorm(1))=(dattot-tot)/totr
	pff=pfsave
c
40	maxit=10
	dellim=5.
	irspec=0
	call manxct(confit)
	call prtpar(.false.)
	fitalp(1)=ftalpo
	ntr=ntrsav
	call sunny
	if (ntr.eq.1) goto 22
	amul=1.+conal(1)-2.*alpha
	param(inorm(2))=param(inorm(1))*(x1a(2)**amul-x2a(2)**amul)/
& 	(x1a(1)**amul-x2a(1)**amul)
c	now let alpha vary
22	if ((.not.fitbmo).and.(.not.ftsigo)) goto 23
	call manxct(confit)
	call prtpar(.false.)
23	fitbem=fitbmo
	fitsga=ftsigo
	if (fitbem) call bmset
	if (fitsga) call sigset
	return
c
c9	para=0.
c	do 14 it=1,notrig
c14	para=para+prbrem(e1(it),ie1(it),ie2(it),de(it),binb(1,it))
c	param(ifixed)=para/float(notrig)
c	return
c
20	iuse=ifix((erup(1)-e1(1))/de(1))+ie1(1)
	call regset(e1(1),erup(1),ptoe(ie1(1),1),ptoe(iuse,1),1)
	iuse=ifix((erdown(2)-e1(1))/de(1))+ie1(1)
	call regset(erdown(2),e2(1),ptoe(iuse,1),ptoe(ie2(1),1),2)
	x1=x1a(1)
	x2=x2a(1)
	goto 21
	end
