Window sampling

window_sampling_sin.nim

## Super dump MC sampling over the entrance window using the Johanna's code from `raytracer2018.nim`
## to check the coverage of the strongback of the 2018 window
import ggplotnim, random, chroma
proc colorMe(y: float): bool =
  const
    stripDistWindow = 2.3  #mm
    stripWidthWindow = 0.5 #mm
  if abs(y) > stripDistWindow / 2.0 and
     abs(y) < stripDistWindow / 2.0 + stripWidthWindow or
     abs(y) > 1.5 * stripDistWindow + stripWidthWindow and
     abs(y) < 1.5 * stripDistWindow + 2.0 * stripWidthWindow:
    result = true
  else:
    result = false

proc sample() =
  randomize(423)
  const nmc = 100_000
  let black = color(0.0, 0.0, 0.0)
  var dataX = newSeqOfCap[float](nmc)
  var dataY = newSeqOfCap[float](nmc)
  var inside = newSeqOfCap[bool](nmc)
  for idx in 0 ..< nmc:
    let x = rand(-7.0 .. 7.0)
    let y = rand(-7.0 .. 7.0)
    if x*x + y*y < 7.0 * 7.0:
      dataX.add x
      dataY.add y
      inside.add colorMe(y)
  let df = seqsToDf(dataX, dataY, inside)
  echo "A fraction of ", df.filter(f{`inside` == true}).len / df.len, " is occluded by the strongback"
  let dfGold = df.filter(f{abs(idx(`dataX`, float)) <= 2.25 and
                           abs(idx(`dataY`, float)) <= 2.25})
  echo "Gold region: A fraction of ", dfGold.filter(f{`inside` == true}).len / dfGold.len, " is occluded by the strongback"
  ggplot(df, aes("dataX", "dataY", fill = "inside")) +
    geom_point() +
    # draw the gold region as a black rectangle
    geom_linerange(aes = aes(y = 0, x = 2.25, yMin = -2.25, yMax = 2.25), color = some(black)) +
    geom_linerange(aes = aes(y = 0, x = -2.25, yMin = -2.25, yMax = 2.25), color = some(black)) +
    geom_linerange(aes = aes(x = 0, y = 2.25, xMin = -2.25, xMax = 2.25), color = some(black)) +
    geom_linerange(aes = aes(x = 0, y = -2.25, xMin = -2.25, xMax = 2.25), color = some(black)) +
    xlab("x [mm]") + ylab("y [mm]") +
    ggsave("/home/basti/org/Figs/statusAndProgress/detector/SiN_window_occlusion.png", width = 1150, height = 1000)
sample()