anchor.diff:  !++  Calculate aperture of collimator
   !JUNE2005   HERE ONE HAS TO HAVE PARTICULAR TREATMENT OF THE OPENING OF
   !JUNE2005   THE PRIMARY COLLIMATOR OF RHIC
  -  if(cdb_cNameUC(icoll)(1:4).ne.'COLM') then
  -    nsig = nsig + gap_rms_error(icoll)
  -    xmax = nsig*sqrt(bx_dist*myemitx0_collgap)
  -    ymax = nsig*sqrt(by_dist*myemity0_collgap)
  -    xmax_pencil = (nsig+pencil_offset)*sqrt(bx_dist*myemitx0_collgap)
  -    ymax_pencil = (nsig+pencil_offset)*sqrt(by_dist*myemity0_collgap)
  -    xmax_nom   = cdb_cNSig(icoll)*sqrt(cdb_cBx(icoll)*myemitx0_collgap)
  -    ymax_nom   = cdb_cNSig(icoll)*sqrt(cdb_cBy(icoll)*myemity0_collgap)
  -    c_rotation = cdb_cRotation(icoll)
  -    c_length   = cdb_cLength(icoll)
  -    c_material = cdb_cMaterial(icoll)
  -    c_offset   = cdb_cOffset(icoll)
  -    c_tilt(1)  = cdb_cTilt(1,icoll)
  -    c_tilt(2)  = cdb_cTilt(2,icoll)
  -
  -    calc_aperture   = sqrt( xmax**2 * cos_mb(c_rotation)**2 + ymax**2 * sin_mb(c_rotation)**2 )
  -    nom_aperture    = sqrt( xmax_nom**2 * cos_mb(c_rotation)**2 + ymax_nom**2 * sin_mb(c_rotation)**2 )
  -    pencil_aperture = sqrt( xmax_pencil**2 * cos_mb(c_rotation)**2+ ymax_pencil**2 * sin_mb(c_rotation)**2 )
  +  nsig = nsig + gap_rms_error(icoll)
  +  xmax = nsig*sqrt(bx_dist*myemitx0_collgap)
  +  ymax = nsig*sqrt(by_dist*myemity0_collgap)
  +  xmax_pencil = (nsig+pencil_offset)*sqrt(bx_dist*myemitx0_collgap)
  +  ymax_pencil = (nsig+pencil_offset)*sqrt(by_dist*myemity0_collgap)
  +  xmax_nom   = cdb_cNSig(icoll)*sqrt(cdb_cBx(icoll)*myemitx0_collgap)
  +  ymax_nom   = cdb_cNSig(icoll)*sqrt(cdb_cBy(icoll)*myemity0_collgap)
  +  c_rotation = cdb_cRotation(icoll)
  +  c_length   = cdb_cLength(icoll)
  +  c_material = cdb_cMaterial(icoll)
  +  c_offset   = cdb_cOffset(icoll)
  +  c_tilt(1)  = cdb_cTilt(1,icoll)
  +  c_tilt(2)  = cdb_cTilt(2,icoll)
  +
  +  calc_aperture   = sqrt( xmax**2 * cos_mb(c_rotation)**2 + ymax**2 * sin_mb(c_rotation)**2 )
  +  nom_aperture    = sqrt( xmax_nom**2 * cos_mb(c_rotation)**2 + ymax_nom**2 * sin_mb(c_rotation)**2 )
  +  pencil_aperture = sqrt( xmax_pencil**2 * cos_mb(c_rotation)**2+ ymax_pencil**2 * sin_mb(c_rotation)**2 )
   
   !++  Get x and y offsets at collimator center point
  -    x_pencil(icoll) = xmax_pencil * (cos_mb(c_rotation))
  -    y_pencil(icoll) = ymax_pencil * (sin_mb(c_rotation))
  +  x_pencil(icoll) = xmax_pencil * (cos_mb(c_rotation))
  +  y_pencil(icoll) = ymax_pencil * (sin_mb(c_rotation))
   
   !++  Get corresponding beam angles (uses xp_max)
  -    xp_pencil(icoll) = -one * sqrt(myemitx0_collgap/tbetax(ie))*talphax(ie)* xmax / sqrt(myemitx0_collgap*tbetax(ie))
  -    yp_pencil(icoll) = -one * sqrt(myemity0_collgap/tbetay(ie))*talphay(ie)* ymax / sqrt(myemity0_collgap*tbetay(ie))
  -    xp_pencil0 = xp_pencil(icoll)
  -    yp_pencil0 = yp_pencil(icoll)
  -
  -    pencil_dx(icoll) = sqrt(xmax_pencil**2 * cos_mb(c_rotation)**2 + ymax_pencil**2 * sin_mb(c_rotation)**2)-calc_aperture
  +  xp_pencil(icoll) = -one * sqrt(myemitx0_collgap/tbetax(ie))*talphax(ie)* xmax / sqrt(myemitx0_collgap*tbetax(ie))
  +  yp_pencil(icoll) = -one * sqrt(myemity0_collgap/tbetay(ie))*talphay(ie)* ymax / sqrt(myemity0_collgap*tbetay(ie))
  +  xp_pencil0 = xp_pencil(icoll)
  +  yp_pencil0 = yp_pencil(icoll)
  +
  +  pencil_dx(icoll) = sqrt(xmax_pencil**2 * cos_mb(c_rotation)**2 + ymax_pencil**2 * sin_mb(c_rotation)**2)-calc_aperture
   
   !++ TW -- tilt for of jaw for pencil beam
   !++ as in Ralphs orig routine, but not in collimate subroutine itself
  
positive.code:     sum      = zero
      sqsum    = zero
  
      do j = 1, napx
        if( part_hit_pos (j).eq.ie .and. part_hit_turn(j).eq.iturn ) then
          if(part_impact(j).lt.-half) then
            write(lout,*) 'ERR>  Found invalid impact parameter!', part_impact(j)
            write(outlun,*) 'ERR>  Invalid impact parameter!', part_impact(j)
            call prror
          end if
  
negative.code: 
  ! should name_sel(1:11) extended to allow longer names as done for
  ! coll the coll_ellipse.dat file !!!!!!!!
    if((cdb_cNameUC(icoll) == name_sel .or. cdb_cName(icoll) == name_sel) .and. iturn == 1) then
      num_selhit = 0
      num_surhit = 0
      num_selabs = 0
  
      do j=1,napx
        if(part_hit_pos (j) == ie .and. part_hit_turn(j) == iturn) then
          num_selhit = num_selhit+1
          if(part_abs_pos(j)  == 0 .and. part_abs_turn(j) == 0) then
            num_surhit = num_surhit+1
          else
            num_selabs = num_selabs + 1
          end if
  !++  If we want to select only partciles interacting at the specified
  !++  collimator then remove all other particles and reset the number
  !++  of the absorbed particles to the selected collimator.
        else if(do_select .and. firstrun) then
          part_select(j) = 0
          n_tot_absorbed = num_selabs
        end if
  