; rate_comparisons_am.pro
;
; test the measured count rates for the block 6 am241 runs
; to see if they match the predicted total count rates

; Do this when printing to x window

red=255L

; Do these when printing to ps file

;r=[0,1]
;g=[0,0]
;b=[0,0]
;set_plot,'ps'
;tvlct,255*r,255*g,255*b
;red=1

; items to change from run to run

filename="blk_6_am241_200_20_0_0_021127_1"
intensity_26=7.96e4   ; am-241-612 26 keV photons/s in 4 pi on 12/01/02
intensity_33=3.72e3   ; am-241-612 33 keV photons/s in 4 pi on 12/01/02
intensity_60=1.14e6   ; am-241-612 60 keV photons/s in 4 pi on 12/01/02
theta=0

; items constant for Am241 runs in block 6

dpi_dir="/local/data/gris1h/block_cal/Block_6/raw/"
title1="am241"
ps_filename=filename+"_rate_comparison.ps"

; calculate source position (in cartesian)

source_x=-100*sin(theta*3.14159/180.)+0.01
source_y=0.01
source_z=100*cos(theta*3.14159/180.)

; print source position

print, "source position:"
print, "  x: ",source_x
print, "  y: ",source_y
print, "  z: ",source_z
print, "  tan(theta): ",sqrt(source_x*source_x+source_y*source_y)/$
      sqrt(source_x*source_x+source_y*source_y+source_z*source_z)

; create a good_map that includes only center detectors

good_map=intarr(34,85)+1
good_map[0,*]=0
good_map[15:18,*]=0
good_map[33,*]=0
for n=0,7 do good_map[*,11*n]=0
for n=0,7 do good_map[*,11*n+7]=0
for n=0,6 do good_map[*,11*n+8]=0
for n=0,6 do good_map[*,11*n+9]=0
for n=0,6 do good_map[*,11*n+10]=0

; make detx and dety (images)

detx=fltarr(34,85)
for i=0,84 do detx[*,i]=findgen(34)
detx=0.42*(temporary(detx)-16.5)

dety=fltarr(34,85)
for i=0,33 do dety[i,*]=findgen(85)
dety=0.42*(temporary(dety)-42)

; make detector flat-fielding correction (cos_factor and r2)

relative_x=source_x-detx
relative_y=source_y-dety
relative_z=source_z

temp1=0.05*abs(relative_z/relative_x)
temp2=0.05*abs(relative_z/relative_y)

cos_factor=$
  cos(atan(sqrt(relative_x*relative_x+relative_y*relative_y)/relative_z))+$
  ((temp1 lt 0.5)*temp1+(temp1 ge 0.5)*0.5)*$
  cos(atan(sqrt(relative_y*relative_y+relative_z*relative_z)/relative_x))+$
  ((temp2 lt 0.5)*temp2+(temp2 ge 0.5)*0.5)*$
  cos(atan(sqrt(relative_x*relative_x+relative_z*relative_z)/relative_y))

r2=relative_x*relative_x+relative_y*relative_y+relative_z*relative_z

; make CZT quantum efficiency factor

top_cos=cos(atan(sqrt(relative_x*relative_x+relative_y*relative_y)/relative_z))
qe_26=1-exp(-54.320*0.2/top_cos)   ; this is pe
;qe_26=1-exp(-54.876*0.2/top_cos)  ; this is pe+incoh (should not be used)
qe_33=1-exp(-173.894*0.2/top_cos)
;qe_33=1-exp(-174.486*0.2/top_cos)
qe_60=1-exp(-36.231*0.2/top_cos)
;qe_60=1-exp(-36.875*0.2/top_cos)

; make attenuation factor

f_26=0.924
f_33=0.948
f_60=0.972
;f_26=1.0
;f_33=1.0
;f_60=1.0

; a couple of other constants

fourpi=12.5663706144
area=0.16
; predicted individual rates

rates_predicted=$
  intensity_26/fourpi/r2*f_26*cos_factor*area*qe_26+$
  intensity_33/fourpi/r2*f_33*cos_factor*area*qe_33+$
  intensity_60/fourpi/r2*f_60*cos_factor*area*qe_60

; measurement

; make arrays of counts and count rates (from dpi)
; note: when making a dpi from a single block calibration file,
; you must first change the BLOCK_ID keywords in the SPECTRUMXXX headers
; to the appropriate number (for block 6, BLOCK_ID = 11)

dpi_big=readfits(dpi_dir+filename+".dpi",exten_no=0)
exposure=(mrdfits(dpi_dir+filename+".dpi",1)).exposure
if (max(exposure) ne min(exposure)) then $
  print, "WARNING: max(exposure) ne min(exposure)" 
print, "min(exposure): ",min(exposure) 
print, "max(exposure): ",max(exposure)
rates_big=dpi_big/max(exposure)
rates=rates_big[108:141,0:84]

; histograms 

;set_plot,'ps'
;device,filename=ps_filename,/color
plot,indgen(400)*0.01,histogram(rates,binsize=0.01,min=0),yrange=[0,500],$
  xtitle="counts/s",$
  title=title1+": Number of detectors per bin!Cbin size: 0.01 counts/s",$
  ymargin=[4,4]
oplot,indgen(400)*0.01,histogram(rates_predicted,binsize=0.01,min=0),color=red
;device,/close
;set_plot,'x'

; construct a "good index table"

good_ind=where((rates*good_map gt 0.5) and (rates*good_map lt 2))

; averages

print,mean(rates[good_ind])
print,mean(rates_predicted[good_ind])

print,"ratio:",mean(rates[good_ind])/mean(rates_predicted[good_ind])

end