* phi_crop.kumac - create a mask for cropping a theta,phi plot
*                  of intensity or polarization, for a given
*                  collimator aperture and beam centroid offset.
*                  If the offset is zero, the mask is just a simple
*                  cutoff in theta, otherwise it is more complicated.
macro make id=90 x0=0 y0=0 coldiam=3.4e-3 dist=75 Emax=12 opt=''
  thetacut=$sigma([coldiam]/2/[dist]*[Emax]/0.511e-3)
  thetax0=$sigma([x0]/[dist]*[Emax]/0.511e-3)
  thetay0=$sigma([y0]/[dist]*[Emax]/0.511e-3)
  fun2 [id] (x*cos(y)-([thetax0]))**2+(x*sin(y)-([thetay0]))**2<[thetacut]**2 100 0 $sigma(2*[thetacut]) 50 0 3.1416 [opt]
return

macro sum id=100 idcrop=90 verbose=1
  if ($hexist([id]).eq.0) then
    message histogram [id] does not exist!
    exitm
  elseif ($hexist([idcrop]).eq.0) then
    exec phi_crop#make
  endif
  idwork=12345
  idwork2=123456
  if ($hexist([idwork]).ne.0) then
    h/del [idwork]
    h/del [idwork2]
  endif
  op/mult [id] [idcrop] [idwork]
  prox [idwork]
  h/proj [idwork]
  h/copy [idwork].prox [idwork2]
  h/del [idwork]
  add [idwork2] [idwork2] [idwork] $sigma(2*pi/50) 0
  v/del rfw*
  v/cre rfw1(100)
  v/cre rfw2(100)
  get/abs [idwork] rfw1
  get/con [idwork] rfw2
  sigma rfw2=rfw2*rfw1
  put/con [idwork] rfw2
  sum=$hinfo([idwork],sum)
  xbins=$hinfo([idwork],xbins)
  xmin=$hinfo([idwork],xmin)
  xmax=$hinfo([idwork],xmax)
  sigma hinteg=[sum]*([xmax]-[xmin])/[xbins]
  if ([verbose].gt.0) then
    message sum is $sigma(hinteg(1))
  endif
return