!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!	ANALYSIS PROGRAM FOR BCM CALIBRATION RUNS 
!
!   (Conversion from Dave J. Mack's PLOTDATA version to PHYSICA)
!
!						author: Heinz C. Anklin
!						date:   may '97
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
! feb-26-98:  - tiled the plotting window so we could get more result
!		plots in a single hclog screen capture
!                                                       djm
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
! july-14-97:  - added several online comments
!              - added plot: residuals versus index
!                                                       Heinz C. Anklin
!
!-------------------------------------------------------------------------
!
DEFAULTS
DESTROY *
DISABLE ECHO
CLEAR ALL
!
runid = `0'
!
gainunse = 250.06E-06
!gainunse = 250.0E-06
gainbcm1 = 280.0E-06		! start value for fit
gainbcm2 = 280.0E-06		! start value for fit
gainbcm3 = 237.0E-06		! start value for fit
!
DISPLAY ` '
DISPLAY `Warning: Be aware that the BCM calibration you have just started'
DISPLAY `         is based on a fixed UNSER monitor gain !'
DISPLAY `         The actual value in units of microA/Hz is: '
=gainunse
DISPLAY `last checked on:' 
DISPLAY `March 99 DJM'
!
DISPLAY ` '
DISPLAY `Have you already copied the file:   charge#####.txt' 
DISPLAY `to current directory. If not, do it now.'
DISPLAY ` '
!
INQUIRE `Enter the runnumber:' runid
offfit = `n'
a = `y'
infile = `charge'//runid//`.txt'
READ infile sample unserate bcm1rate bcm2rate bcm3rate interval
!
nmax = LEN(sample)
bindur= NINT(interval[2])
!
!******************************  Plot Rates   **************************
!
!
SET
 cursor -2
 %xloc 50
 %yloc 92.5
 pchar 0
 charsz 0.50
 color 1
 txthit 0.7

AUTOSCALE ON
LABEL\xaxis `Sample'
LABEL\yaxis `Rate (Hz)'
TEXT `Run:  '//runid 
Set txthit 0.70
LEGEND ON
LEGEND FRAME OFF
LEGEND AUTOHEIGHT OFF
LEGEND frame 15 73 25 90
GRAPH `Unser' sample unserate
SET color 2
GRAPH\noaxes `BCM1' sample bcm1rate
SET color 3
GRAPH\noaxes `BCM2' sample bcm2rate
SET color 5
GRAPH\noaxes `BCM3' sample bcm3rate
CLEAR\noreplot
REPLOT
LEGEND OFF
!
DISPLAY ` '
DISPLAY `Plot 1: Rates of the UNSER monitor BCM1, BCM2 and BCM3 versus sample number.'
DISPLAY `        Each sample corresponds to an interval of [s]:'
=bindur
DISPLAY ` '
!
!
!********************   Average Beam OFF intervals  *********************
!
BEAM_OFF:
!
DISPLAY ` '
DISPLAY `Enter pairs of points delineating BEAM OFF intervals.'
DISPLAY `Begin at small sample numbers and work chronically to the end'
DISPLAY `of the run.'
DISPLAY ` '
DISPLAY `Make sure that no beam spikes are included in the'
DISPLAY `selected intervals.'
DISPLAY `Press Q when finished.'
DISPLAY ` '
!
PICK xu yu
ixu = int(xu)
L = LEN(ixu)
e1 = L/2
e2 = INT(L/2)
IF (e1~=e2) THEN
  DISPLAY `Even number of points needed!'
  DISPLAY `TRY AGAIN!!!'
  DESTROY xu yu
  TERMINAL
  GOTO BEAM_OFF
ENDIF
!			
ic = 0
DO i = [1:L:2]
 ic = ic + 1
 nbin[ic] = ixu[i+1] - ixu[i]
 STATISTICS unserate[ixu[i]:ixu[i+1]] unsezero[ic]\mean unsezstd[ic]\sdev
 STATISTICS bcm1rate[ixu[i]:ixu[i+1]] bcm1zero[ic]\mean bcm1zstd[ic]\sdev
 STATISTICS bcm2rate[ixu[i]:ixu[i+1]] bcm2zero[ic]\mean bcm2zstd[ic]\sdev
 STATISTICS bcm3rate[ixu[i]:ixu[i+1]] bcm3zero[ic]\mean bcm3zstd[ic]\sdev
ENDDO 
!
! average all specified beam off intervals:
!
SCALAR\DUMMY J
totbin[1] = SUM(nbin[J],J,1:e1)
ndur = nbin * bindur
Zunser    = SUM(unsezero[J]*nbin[J],J,1:e1) / totbin[1]
bcm1[1]   = SUM(bcm1zero[J]*nbin[J],J,1:e1) / totbin[1]
bcm2[1]   = SUM(bcm2zero[J]*nbin[J],J,1:e1) / totbin[1]
bcm3[1]   = SUM(bcm3zero[J]*nbin[J],J,1:e1) / totbin[1]
!
offsetbcm1 = bcm1[1]
offsetbcm2 = bcm2[1]
offsetbcm3 = bcm3[1]
!
!***************  Average Beam ON intervals **************************
!
BEAM_ON:
!
DISPLAY ` '
DISPLAY `Enter pairs of points delineating BEAM ON intervals.'
DISPLAY `Begin at small sample numbers and work chronically to the end'
DISPLAY `of the run.'
DISPLAY ` '
DISPLAY `Make sure that only stable beam periods are included in the'
DISPLAY `selected intervals.'
DISPLAY `Press Q when finished.'
DISPLAY ` '
!
PICK xu yu
ixu = int(xu)
L = LEN(ixu)
!
e1 = L/2
e2 = INT(L/2)
IF (e1~=e2) THEN
  DISPLAY `Even number of points needed!'
  DISPLAY `TRY AGAIN!!!'
  DESTROY xu yu
  TERMINAL
  GOTO BEAM_ON
ENDIF
!
Iunser[1]  = 0.0
is = 1
!
unsestd[1] = 0.0
DO i = [1:L:2]
 is = is + 1
 STATISTICS unserate[ixu[i]:ixu[i+1]] unser\mean unsestd[is]\sdev
 STATISTICS unserate[ixu[i]:ixu[i+1]] unservec[is]\mean unsestd[is]\sdev
 STATISTICS bcm1rate[ixu[i]:ixu[i+1]] bcm1[is]\mean bcm1std[is]\sdev
 STATISTICS bcm2rate[ixu[i]:ixu[i+1]] bcm2[is]\mean bcm2std[is]\sdev
 STATISTICS bcm3rate[ixu[i]:ixu[i+1]] bcm3[is]\mean bcm3std[is]\sdev
 Iunser[is] = gainunse * (unser - Zunser)
 totbin[is] = ixu[i+1] - ixu[i]
 indx[is] = is - 1
ENDDO
!
intdur = totbin * bindur
!
!
!******************  OUTPUT to file  ***********************************
!
fileout = `run'//runid//`.cal'
WRITE\SCALAR\FORMAT fileout  ('run:', f8.0) runid
WRITE\TEX\APPEND fileout `unser (Hz)'
WRITE\APPEND\FORMAT fileout (F13.1)  unsezero
WRITE\APPEND\FORMAT fileout (F13.1)  unservec

WRITE\TEX\APPEND fileout `     Current(muA)  BCM1(Hz) BCM2(Hz)  BCM3(Hz)     Weight'
WRITE\APPEND\FORMAT fileout (6F13.1)  iunser  bcm1 bcm2 bcm3 totbin
!
!
INPUT2:
!
answ  = `y'
INQUIRE `REMOVE I = 0 from fit (RECOMMENDED!)? <y>/n:' answ
IF EQS(answ,`n') THEN GOTO JUMP
IF EQS(answ,`N') THEN GOTO JUMP
IF NES(answ,`y') THEN
  IF NES(answ,`Y') THEN
    DISPLAY `invalid key!'
    GOTO INPUT2
  ENDIF
ENDIF  
DESTROY Iunser[1] bcm1[1] bcm2[1] bcm3[1] totbin[1] intdur[1] unsestd[1]
!
JUMP:
!
offfit = `y'
INQUIRE `Fit offset (RECOMMENDED)? <y>/n' offfit
!
! note in outputfile:
!
IF EQS(offfit,`n') THEN WRITE\TEXT\APPEND outfile `** offset not fitted! **'
IF EQS(offfit,`N') THEN WRITE\TEXT\APPEND outfile `** offset not fitted! **'
!
!
!**********************  Graph  UNSER zero intervals *******************
!
UNSER_ZERO:
!
CLEAR
SET 
 color 1
 pchar -3

!UNSER ZERO GETS A WHOLE PAGE TO ITSELF
WINDOW 0
!WINDOW 5
SET txthit 0.7
TEXT `Unser zero Run:'//runid 
LABEL\xaxis `selected interval'
LABEL\yaxis `Residual [micro A]'
GENERATE x 1 1 LEN(unsezero)
yz = (unsezero-Zunser)*gainunse
yzerr = gainunse*unsezstd/sqrt(nbin)
GRAPH  x yz yzerr
ZEROLINES\HORIZONTAL
DISPLAY ` '
DISPLAY `Residuals of the Unser zeroes versus order of the previously  '
DISPLAY `selected intervals. The error bars are statistical only. A    '
DISPLAY `random variation of .2 microA (sigma) is expected. Check for   '
DISPLAY `systematic drifts!'
DISPLAY `Note: The UNSER offset is kept constant at [Hz](wgt. aver.):'
=Zunser
TERMINAL
!
!
!
!***********************  calculate ERRORS  *****************************

!error  = sqrt(1/intdur)	! assuming 1 microA per one second interval
error1 = unsestd/sqrt(totbin) * gainunse     ! standarddeviation from data
!
!
!***************** Fit the data and calculate residuals **********************
!
! BCM1 
!
SCALAR\VARY gainbcm1
SCALAR\VARY offsetbcm1
IF EQS(offfit,`n') THEN SCALAR offsetbcm1
IF EQS(offfit,`N') THEN SCALAR offsetbcm1
FIT\WEIGHTS totbin Iunser = gainbcm1 * (bcm1 - offsetbcm1)
FIT\UPDATE fit1
! Calculate residuals
res1 = fit1 - Iunser
!
! BCM2   
!
SCALAR\VARY gainbcm2
SCALAR\VARY offsetbcm2
IF EQS(offfit,`n') THEN SCALAR offsetbcm2
IF EQS(offfit,`N') THEN SCALAR offsetbcm2
FIT\WEIGHTS totbin Iunser = gainbcm2 * (bcm2 - offsetbcm2)
FIT\UPDATE fit2
! Calculate residuals
res2 = fit2 - Iunser
!
!
! BCM3  
!
SCALAR\VARY gainbcm3
SCALAR\VARY offsetbcm3
IF EQS(offfit,`n') THEN SCALAR offsetbcm3
IF EQS(offfit,`N') THEN SCALAR offsetbcm3
FIT\WEIGHTS totbin Iunser = gainbcm3 * (bcm3 - offsetbcm3)
FIT\UPDATE fit3
! Calculate residuals
res3 = fit3 - Iunser
!
! ********************* Graph residuals versus current ********************
DISPLAY ` '
DISPLAY `Linear fits of BCM rates versus UNSER current have been done.'
DISPLAY `The next plot shows the residuals for BCM1, BCM2, AND BCM3 .'
DISPLAY ` '

CLEAR
WINDOW 1
SET 
 color 1
 pchar -2

LABEL\xaxis `Current (microA)'
LABEL\yaxis `Residual (microA)'
TEXT `Residuals  Run:  '//runid
SET color 2
SET pchar -1
GRAPH Iunser  res1  error1
ZEROLINES\HORIZONTAL
!
!!CLEAR
!WINDOW 6
SET color 3
SET pchar -2
!LABEL\xaxis `Current (microA)'
!LABEL\yaxis `Residual (microA)'
!TEXT `BCM2 Run:  '//runid
GRAPH\noaxes Iunser  res2  error1
!
!!CLEAR
!WINDOW 8
SET color 5
SET pchar -3
!LABEL\xaxis `Current (microA)'
!LABEL\yaxis `Residual (microA)'
!TEXT `BCM3 Run:  '//runid
GRAPH\noaxes Iunser res3  error1
ZEROLINES\HORIZONTAL
!
!
! *********************   Graph Residuals versus time ********************
!

SET
 color 1
 txthit 0.70
 pchar -1

WINDOW 2
!WINDOW 5
LABEL\xaxis `Interval'
LABEL\yaxis `Residual (microA)'
TEXT `Residuals  Run:  '//runid 
Set txthit 0.35
LEGEND ON
LEGEND FRAME OFF
LEGEND AUTOHEIGHT OFF
LEGEND frame 15 73 25 90
!GRAPH `Unser' sample unsercur
SET color 2
SET pchar -1
GRAPH        `BCM1' indx res1 error1
ZEROLINES\HORIZONTAL
SET color 3
SET pchar -2
GRAPH\noaxes `BCM2' indx res2 error1
SET color 5
SET pchar -3
GRAPH\noaxes `BCM3' indx res3 error1
!CLEAR\noreplot
!REPLOT
LEGEND OFF
!
DISPLAY ` '
DISPLAY `Check for systematic drifts in the residuals versus the'
DISPLAY `order of the selected intervals!'
DISPLAY ` '
TERMINAL
!
!
!**********************    GRAPH RATES AGAIN   *************************
CLEAR
SET
 cursor -2
 %xloc 50
 %yloc 92.5
 pchar 0
 charsz 0.50
 color 1
 txthit 0.70

AUTOSCALE ON
WINDOW 1
LABEL\xaxis `Sample'
LABEL\yaxis `Rate (Hz)'
TEXT `Run:  '//runid 
Set txthit 0.70
LEGEND ON
LEGEND FRAME OFF
LEGEND AUTOHEIGHT OFF
LEGEND frame 15 73 25 90
GRAPH `Unser' sample unserate
SET color 2
GRAPH\noaxes `BCM1' sample bcm1rate
SET color 3
GRAPH\noaxes `BCM2' sample bcm2rate
SET color 5
GRAPH\noaxes `BCM3' sample bcm3rate
!?? CLEAR\noreplot
REPLOT
LEGEND OFF

!**********************    Graph all currents   *************************
!
unsercur = gainunse * (unserate - Zunser)
bcm1cur = gainbcm1 * (bcm1rate - offsetbcm1)
bcm2cur = gainbcm2 * (bcm2rate - offsetbcm2)
bcm3cur = gainbcm3 * (bcm3rate - offsetbcm3)
SET
 color 1
 txthit 0.70
 pchar 0

!CLEAR
WINDOW 2
LABEL\xaxis `Sample'
LABEL\yaxis `Current (microA)'
TEXT `Run:  '//runid 
Set txthit 0.70
LEGEND ON
LEGEND FRAME OFF
LEGEND AUTOHEIGHT OFF
LEGEND frame 15 73 25 90
GRAPH `Unser' sample unsercur
SET color 2
GRAPH\noaxes `BCM1' sample bcm1cur
SET color 3
GRAPH\noaxes `BCM2' sample bcm2cur
SET color 5
GRAPH\noaxes `BCM3' sample bcm3cur
!CLEAR\noreplot
!REPLOT
LEGEND OFF
!
DISPLAY ` '
DISPLAY `Rates and current from the calibration run for BCM1,2,3. Due to     '
DISPLAY `nonlinearities at very low currents a large discrepancy is expected  '
DISPLAY `for beam off intervals (particularly BCM3).'
DISPLAY ` '
!
TERMINAL
!
!
DISPLAY ` '
INQUIRE `See the plots again? <y>/n' a
IF EQS(a,`y') THEN GOTO UNSER_ZERO
IF EQS(a,`Y') THEN GOTO UNSER_ZERO
!
!
!******************************  OUTPUT  ******************************
!
!
WRITE\TEXT\APPEND fileout `gain factors for three cavity monitors'
WRITE\SCALAR\APPEND\FORMAT fileout ('gbcm1_gain = ',f10.8,-
'              ; microA/Hz') gainbcm1
!
WRITE\SCALAR\APPEND\FORMAT fileout ('gbcm2_gain = ',f10.8,-
'              ; microA/Hz') gainbcm2
!
WRITE\SCALAR\APPEND\FORMAT fileout ('gbcm3_gain = ',f10.8,-
'              ; microA/Hz') gainbcm3
!
WRITE\SCALAR\APPEND\FORMAT fileout ('gunser_gain = ',f10.8,-
'              ; microA/Hz') gainunse

WRITE\TEXT\APPEND fileout `zero offsets for BCM s'
!
WRITE\SCALAR\APPEND\FORMAT fileout ('gbcm1_offset = ',f8.0,-
'              ; Hz') offsetbcm1
!
WRITE\SCALAR\APPEND\FORMAT fileout ('gbcm2_offset = ',f8.0,-
'              ; Hz') offsetbcm2
!
WRITE\SCALAR\APPEND\FORMAT fileout ('gbcm3_offset = ',f8.0,-
'              ; Hz') offsetbcm3
!
WRITE\SCALAR\APPEND\FORMAT fileout ('gunser_offset = ',f8.0,-
'              ; Hz') Zunser
!
DISPLAY ` '
DISPLAY `The fitted gains (and offsets) have been written to file:'
=fileout
DISPLAY ` '
!
!WRITE\SCALAR\APPEND\FORMAT fileout ('BCM1  gain[nA/kHz]: ',f10.1,-
!'  offset[kHz]: ',f12.3) gainbcm1*1e06 offsetbcm1/1E03
!
!WRITE\SCALAR\APPEND\FORMAT fileout ('BCM2  gain[nA/kHz]: ',f10.1,-
!'  offset[kHz]: ',f12.3) gainbcm2*1E06 offsetbcm2/1E03
!
!WRITE\SCALAR\APPEND\FORMAT fileout ('BCM3  gain[nA/kHz]: ',f10.1,-
!'  offset[kHz]: ',f12.3) gainbcm3*1E06 offsetbcm3/1E03