anchor.diff:              by_dist = tbetay(j)
             end if
   
  -          sig_offset = cdb_cOffset(i)/(sqrt(bx_dist**2 * cos_mb(cdb_cRotation(i))**2 + by_dist**2 * sin_mb(cdb_cRotation(i))**2 ))
  -          write(coll_twissLikeUnit,*) bez(myix),tbetax(j),tbetay(j), torbx(j),torby(j), nsig, gap_rms_error(i)
  -          write(coll_sigmaSetUnit,*) bez(myix), gap_h1, gap_h2, gap_h3, gap_h4, sig_offset, cdb_cOffset(i), nsig, gap_rms_error(i)
  -
  -          if((gap_h1 + sig_offset) .le. mingap) then
  -            mingap = gap_h1 + sig_offset
  +          sig_offset = cdb_cOffset(i)/(sqrt(bx_dist**2 * cos_mb(cdb_cRotation(i))**2 + by_dist**2 * sin_mb(cdb_cRotation(i))**2))
  +          write(coll_twissLikeUnit,*) cdb_cName(i),tbetax(j),tbetay(j),torbx(j),torby(j),nsig,gap_rms_error(i)
  +          write(coll_sigmaSetUnit,*) cdb_cName(i),gap_h1,gap_h2,gap_h3,gap_h4,sig_offset,cdb_cOffset(i),nsig,gap_rms_error(i)
  +
  +          if((gap_h1 + sig_offset) <= mingap) then
  +            mingap         = gap_h1 + sig_offset
               coll_mingap_id = i
  -            coll_mingap2 = cdb_cName(i)
  -          else if((gap_h2 + sig_offset) .le. mingap) then
  -            mingap = gap_h2 + sig_offset
  +            coll_mingap2   = cdb_cName(i)
  +          else if((gap_h2 + sig_offset) <= mingap) then
  +            mingap         = gap_h2 + sig_offset
               coll_mingap_id = i
  -            coll_mingap2 = cdb_cName(i)
  -          else if((gap_h3 - sig_offset) .le. mingap) then
  -            mingap = gap_h3 - sig_offset
  +            coll_mingap2   = cdb_cName(i)
  +          else if((gap_h3 - sig_offset) <= mingap) then
  +            mingap         = gap_h3 - sig_offset
               coll_mingap_id = i
  -            coll_mingap2 = cdb_cName(i)
  -          else if((gap_h4 - sig_offset) .le. mingap) then
  -            mingap = gap_h4 - sig_offset
  +            coll_mingap2   = cdb_cName(i)
  +          else if((gap_h4 - sig_offset) <= mingap) then
  +            mingap         = gap_h4 - sig_offset
               coll_mingap_id = i
  -            coll_mingap2 = cdb_cName(i)
  +            coll_mingap2   = cdb_cName(i)
             end if
           end if
         end if
  
positive.code:     end if
  
      do i=1,cdb_nColl
  ! start searching minimum gap
        if(cdb_cName(i) == bez(myix)) then
          if( cdb_cLength(i) > zero ) then
            nsig_err = nsig + gap_rms_error(i)
  
  ! jaw 1 on positive side x-axis
            gap_h1 = nsig_err - sin_mb(cdb_cTilt(1,i))*cdb_cLength(i)/2
            gap_h2 = nsig_err + sin_mb(cdb_cTilt(1,i))*cdb_cLength(i)/2
  
  ! jaw 2 on negative side of x-axis (see change of sign comapred
  ! to above code lines, alos have a look to setting of tilt angle)
            gap_h3 = nsig_err + sin_mb(cdb_cTilt(2,i))*cdb_cLength(i)/2
            gap_h4 = nsig_err - sin_mb(cdb_cTilt(2,i))*cdb_cLength(i)/2
  
  ! find minumum halfgap
  ! --- searching for smallest halfgap
  !! ---scaling for beta beat needed?
  !                        if (do_nominal) then
  !                           bx_dist = cdb_cBx(icoll) * scale_bx / scale_bx0
  !                           by_dist = cdb_cBy(icoll) * scale_by / scale_by0
  !                        else
  !                           bx_dist = tbetax(j) * scale_bx / scale_bx0
  !                           by_dist = tbetay(j) * scale_by / scale_by0
  !                        endif
            if (do_nominal) then
              bx_dist = cdb_cBx(icoll)
              by_dist = cdb_cBy(icoll)
            else
  
negative.code:   rnd_k2  = 0
    call rluxgo(rnd_lux, c_offsettilt_seed, rnd_k1, rnd_k2)
  
    ! Generate random tilts (Gaussian distribution plus systematic)
    if(c_rmstilt_prim > zero .or. c_rmstilt_sec > zero .or. c_systilt_prim /= zero .or. c_systilt_sec /= zero) then
      do icoll = 1, cdb_nColl
        if(cdb_cType(icoll) == cdb_typPrimary) then
          c_rmstilt = c_rmstilt_prim
          c_systilt = c_systilt_prim
        else
          c_rmstilt = c_rmstilt_sec
          c_systilt = c_systilt_sec
        end if
        cdb_cTilt(1,icoll) = c_systilt+c_rmstilt*ran_gauss2(three)
        if(systilt_antisymm) then
          cdb_cTilt(2,icoll) = -one*c_systilt+c_rmstilt*ran_gauss2(three)
        else
          cdb_cTilt(2,icoll) =      c_systilt+c_rmstilt*ran_gauss2(three)
        end if
        write(outlun,*) 'INFO>  Collimator ',trim(cdb_cName(icoll)),' jaw 1 has tilt [rad]: ',cdb_cTilt(1,icoll)
        write(outlun,*) 'INFO>  Collimator ',trim(cdb_cName(icoll)),' jaw 2 has tilt [rad]: ',cdb_cTilt(2,icoll)
      end do
  
      do icoll=1,cdb_nColl
        if(cdb_cType(icoll) == cdb_typPrimary) then
          cdb_cOffset(icoll) = c_sysoffset_prim + c_rmsoffset_prim*ran_gauss2(three)
        else
          cdb_cOffset(icoll) = c_sysoffset_sec +  c_rmsoffset_sec*ran_gauss2(three)
        end if
        write(outlun,*) 'INFO>  Offset: ',trim(cdb_cName(icoll)),cdb_cOffset(icoll)
      end do
    end if
    do icoll=1,cdb_nColl
      gap_rms_error(icoll) = c_rmserror_gap * ran_gauss2(three)
      write(outlun,*) 'INFO>  Gap RMS error: ',trim(cdb_cName(icoll)),gap_rms_error(icoll)
    end do
  
    ! In case we're using old type jaw fit, this is where we generate the parameters for the new method
    ! After this, the number of slices is also stored per collimator, and can be extracted again later
    call cdb_setMasterJawFit(n_slices, smin_slices, smax_slices, recenter1, recenter2, jaw_fit, jaw_ssf)
  
  !---- creating a file with beta-functions at TCP/TCS
    mingap = 20
  
    do j=1,iu
      ! this transformation gives the right marker/name to the corresponding beta-functions or vice versa
      if(ic(j) <= nblo) then
        myix = mtyp(ic(j),mel(ic(j)))
      else
  