vec/del *
type=all
* set getpoints= 1 if one wants to set the lower and upper gates to define beam current periods.
getpoints=0
if ([type] .eq. '72317') then
npts=7
vec/create fcup([npts]) r 85. 100.  120. 140. 160. 180. 200.
vec/read evid,time,bcm1,bcm2,bcm1hp,bcm2hp,timehp,bcm1hm,bcm2hm,timehm,s503  ../scalers/72317.scal
vec/create gatelo([npts]) r   10.8453  94.5087 178.172 255.638 339.302 453.951 568.601
vec/create gatehi([npts]) r   23.2398  113.101 206.06  277.329 373.387 491.135 636.771
endif
if ([type] .eq. 'all') then
npts=9
vec/create fcup([npts]) r 50. 65. 85. 100.  120. 140. 160. 180. 200.
vec/read evid,time,bcm1,bcm2,bcm1hp,bcm2hp,timehp,bcm1hm,bcm2hm,timehm,s503  ../scalers/all_bcm_calib.scal
vec/create gatelo([npts]) r   37.245  126.716  189.138  272.367  353.515  430.502  515.812  632.333  734.289
vec/create gatehi([npts]) r   76.779  137.120  209.945  299.417  393.049  465.875  563.669  682.270  836.244
endif
sigma stime=sumv(time)
len=$vlen(stime,1)
* online parameters
sigma cur1=(bcm1-250580.)*.01278
sigma cur2=(bcm2-250500.)*.01238
*
*
zone 1 2
graph [len] stime cur1
graph [len] stime cur2
wait
zone 1 1
vec/creat ratbcm1([npts]) r
vec/creat ratbcm2([npts]) r
vec/creat rbcmm2([npts]) r
*
*
zone 1 1
null 0 1400 250000 270000
*graph [len] stime bcm2
pline [len] stime bcm2
if ( [getpoints] .eq. 1) then
vec/create gatelo([npts]) r  
vec/create gatehi([npts]) r  
do i=1,[npts]
fc=fcup([i])
message set low and high gate for [fc] current
vlocate x y s
vec/input gatelo([i]) x(1)
vec/input gatehi([i]) x(2)
enddo
mess time cuts
vec/write gatelo,gatehi ' ' [npts](1x,f8.2)
endif
do i=1,[npts]
lo=gatelo([i])
hi=gatehi([i])
*dline [lo] [lo]  250000 270000 
*dline [hi] [hi]  250000 270000 
enddo
atitle 'Time (sec)' ' Rate ' ! 220
*
do nc=1,[npts]
sb1 = 0
sb2 = 0
ns = 0
 xlo=gatelo([nc])
 xhi=gatehi([nc])
message looping
do i=1,[len]
 xt=stime([i])
 if ( [xt] .lt. [xhi] ) then
	if ( [xt] .gt. [xlo]) then
	  sb1 = [sb1] + bcm1([i])
	  sb2 = [sb2] + bcm2([i])
         ns = [ns]+ 1
        endif
 endif
enddo
fc=fcup([nc])
avebcm1=[sb1]/[ns]
avebcm2=[sb2]/[ns]
vec/input ratbcm1([nc]) [avebcm1]
vec/input ratbcm2([nc]) [avebcm2]
message [nc] [xlo] [xhi]  fc = [fc] bcm1 rate = [avebcm1] bcm2 rate = [avebcm2]  
SET PLCI 2
dline [xlo] [xhi]  [avebcm2] [avebcm2] 
enddo 
SET PLCI 1
vec/write fcup,ratbcm1,ratbcm2 ' ' 3(1x,g15.5)
wait
set mscf 1.0
set mtyp 2.0
vec/cr b1par(2) r 250000. .015
zone 1 2
null 45. 205. 250000 270000
SET MTYP 3
SET PMCI 2
pmarker [npts] fcup ratbcm1 
atitle 'Current (nA)' ' BCM1 rate' ! 220
vec/fit  fcup ratbcm1 ? p1 s 2 b1par
vec/write fcup,ratbcm1 ' ' 2(1x,f15.5)
scal1=1./b1par(2)
con1=b1par(1)
sigma calbcm1=[scal1]*(ratbcm1-[con1])
sigma dbcm1=calbcm1-fcup
message bcm1 constant = [con1] scale = [scal1]
*
vec/cr b2par(2) r 250000. .015
null 45. 205. 250000 270000
pmarker [npts] fcup ratbcm2 apw
atitle 'Current (nA)' ' BCM1 rate' ! 220
vec/fit  fcup ratbcm2 ? p1 s 2 b2par
vec/write fcup,ratbcm2 ' ' 2(1x,f15.5)
scal2=1./b2par(2)
con2=b2par(1)
message bcm2 constant = [con2] scale = [scal2]
wait
sigma calbcm2=[scal2]*(ratbcm2-[con2])
sigma dbcm2=calbcm2-fcup
null 45. 205. -2. 2.
pmarker [npts] fcup dbcm1 apw
atitle 'Current (nA)' ' BCM1 Current - Faraday Cup' ! 220
null 45. 205. -2. 2.
pmarker [npts] fcup dbcm2 apw
atitle 'Current (nA)' ' BCM2 Current - Faraday Cup' ! 220
wait
*sigma dnew1=(ratbcm1-[con1])*[scal1]/((ratbcm1-250580.)*.01278)
*sigma dnew2=(ratbcm2-[con2])*[scal2]/((ratbcm2-250500.)*.01238)
sigma dnew1=(ratbcm1-[con1])*[scal1]/((ratbcm1-252580.)*[scal1])
sigma dnew2=(ratbcm2-[con2])*[scal2]/((ratbcm2-252580.)*[scal2])
graph [npts] fcup dnew1 apw
graph [npts] fcup dnew2 apw
wait
*