! !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! 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 ! !------------------------------------------------------------------------- ! enable echo defaults destroy * clear all device\colour postscript A set font roman.3 set nsxinc 2 set nsyinc 5 set linthk 3 !disable echo ! runid = `0' ! gainunse = 2.498E-04 !2004 !gainunse = 250.01E-06 !2003 !gainunse = 250.06E-06 !summer 99 !gainunse = 250.06E-06 !spring 99 !gainunse = 250.0E-06 !nominal value 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 `JUNE 04 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 #1 ******************************* !****************************** Rates ******************************* !*********************************************************************** ! window 0 ! 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 `Scaler Rates Run: '//runid Set txthit 0.70 ! LEGEND ON LEGEND FRAME OFF LEGEND AUTOHEIGHT OFF LEGEND frame 15 70 35 90 ! GRAPH `Unser' sample unserate ! SET color 2 GRAPH\noaxes `BCM1' sample bcm1rate ! SET color 3 GRAPH\noaxes `BCM2' sample bcm2rate ! CLEAR\noreplot REPLOT terminal ! !******************** Digitize and Average the Beam OFF intervals ********************* ! DISPLAY ` ' DISPLAY `Plot 1: Rates of the UNSER monitor, BCM1, and BCM2 versus sample number.' DISPLAY ` Each sample corresponds to an interval of [s]:' =bindur DISPLAY ` ' ! ! 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] ! !*************** Digitize and Average the 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 digitized means to a 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) Weight' WRITE\APPEND\FORMAT fileout (6F13.1) iunser bcm1 bcm2 totbin ! ! INPUT2: ! answ = `y' INQUIRE `REMOVE I = 0 from fit (RECOMMENDED!)? /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] totbin[1] intdur[1] unsestd[1] ! JUMP: ! offfit = `y' INQUIRE `Fit offset (RECOMMENDED)? /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! **' ! !************************************************************************* !********************** Plot #2 ******************************* !******************* UNSER zero residuals ******************************** !************************************************************************* ! window 0 ! UNSER_ZERO: ! CLEAR set color 1 txthit 0.7 TEXT `Unser Residuals Run:'//runid LABEL\xaxis `Digitized Interval' LABEL\yaxis `Residual (A)' ! GENERATE x 1 1 LEN(unsezero) yz = (unsezero-Zunser)*gainunse yzerr = 2./sqrt(intdur) ! assuming 2 microA per one second interval !yzerr = gainunse*unsezstd/sqrt(nbin) scale 0 20 10 -1 1 4 legend off set color 1 set pchar -8 graph\axes ZEROLINES\HORIZONTAL ! set color 2 GRAPH\noaxes x yz yzerr ! ! 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 = 2./sqrt(intdur) ! assuming 2 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 ! !******************************************************************** !***************************** Plot #3 **************************** !********************* bcm 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.' DISPLAY ` ' CLEAR window 1 SET color 1 pchar -8 LABEL\xaxis `Beam Current (A)' LABEL\yaxis `BCM Fit Residuals (A)' TEXT `BCM Residuals' ! scale 0 100 5 -1 1 4 legend on LEGEND frame 20 15 40 35 graph\axes ZEROLINES\HORIZONTAL ! set color 2 GRAPH\noaxes `BCM1' Iunser res1 error ! SET color 3 GRAPH\noaxes `BCM2' Iunser res2 error ! ! !******************************************************************* !************************** Plot #4 ******************************* !******************** bcm residuals versus time ******************** !******************************************************************* ! window 2 SET color 1 pchar -8 LABEL\xaxis `Digitized Interval' !LABEL\yaxis `BCM Fit Residuals (A)' LABEL\yaxis TEXT `Run: '//runid ! LEGEND ON LEGEND AUTOHEIGHT OFF LEGEND frame 20 15 40 35 ! scale 0 20 4 -1 1 4 !GRAPH `Unser' sample unsercur SET color 1 graph\axes ZEROLINES\HORIZONTAL ! set color 2 GRAPH\noaxes `BCM1' indx res1 error ! SET color 3 GRAPH\noaxes `BCM2' indx res2 error ! DISPLAY ` ' DISPLAY `Check for systematic drifts in the residuals versus the' DISPLAY `order of the selected intervals!' DISPLAY ` ' ! TERMINAL ! !********************************************************************* !************************** Plot #5 ********************************* !********************** GRAPH RATES AGAIN ************************* !********************************************************************* ! CLEAR window 0 SET cursor -2 %xloc 50 %yloc 92.5 pchar 0 charsz 0.50 color 1 txthit 0.70 scale 0 3000 6 0 1000000 5 ! LABEL\xaxis `Sample' LABEL\yaxis `Raw Rate (Hz)' TEXT `Run: '//runid Set txthit 0.70 ! LEGEND ON LEGEND FRAME OFF LEGEND AUTOHEIGHT OFF LEGEND frame 15 70 35 90 ! GRAPH `Unser' sample unserate ! SET color 2 GRAPH\noaxes `BCM1' sample bcm1rate ! SET color 3 GRAPH\noaxes `BCM2' sample bcm2rate ! LEGEND OFF terminal ! !************************************************************************ !************************** Plot #6 and #7 ****************************** !********************* Graph bcm1 and bcm2 currents **************************** !************************************************************************ ! window 0 clear ! unsercur = gainunse * (unserate - Zunser) bcm1cur = gainbcm1 * (bcm1rate - offsetbcm1) bcm2cur = gainbcm2 * (bcm2rate - offsetbcm2) !bcm3cur = gainbcm3 * (bcm3rate - offsetbcm3) SET color 1 txthit 0.70 pchar 0 LABEL\xaxis `Sample' LABEL\yaxis `Current (A)' TEXT `Calibrated Currents' Set txthit 0.50 ! LEGEND ON LEGEND FRAME OFF LEGEND AUTOHEIGHT OFF LEGEND frame 15 70 35 90 ! scale 0 3000 6 -10 100 11 set color 1 GRAPH `Unser' sample unsercur graph\axes ! SET color 2 GRAPH\noaxes `BCM1' sample bcm1cur ! SET color 3 GRAPH\noaxes `BCM2' sample bcm2cur ! DISPLAY ` ' DISPLAY `Rates and current from the calibration run for BCM1. Due to ' DISPLAY `nonlinearities at very low currents a large discrepancy is expected ' DISPLAY `for beam off intervals.' DISPLAY ` ' ! TERMINAL ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! DISPLAY ` ' INQUIRE `See the plots again? /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