2.1.1. Outdated numbers

./writeRunList -b ~/CastData/data/DataRuns2017_Reco.h5 -c ~/CastData/data/CalibrationRuns2017_Reco.h5

./writeRunList -b ~/CastData/data/DataRuns2018_Reco.h5 -c ~/CastData/data/CalibrationRuns2018_Reco.h5

2.7.1. Run 189

The below is copied from thesis.org.

I copied the X-ray finger runs from tpc19 over to ./../../CastData/data/XrayFingerRuns/. The run of interest is mainly the run 189, as it's the run done with the detector installed as in 2017/18 data taking.

cd /dev/shm # store here for fast access & temporary
cp ~/CastData/data/XrayFingerRuns/XrayFingerRun2018.tar.gz .
tar xzf XrayFingerRun2018.tar.gz
raw_data_manipulation -p Run_189_180420-09-53 --runType xray --out xray_raw_run189.h5
reconstruction -i xray_raw_run189.h5 --out xray_reco_run189.h5 # make sure `config.toml` for reconstruction uses `default` clustering!
reconstruction -i xray_reco_run189.h5 --only_charge
reconstruction -i xray_reco_run189.h5 --only_gas_gain
reconstruction -i xray_reco_run189.h5 --only_energy_from_e
plotData --h5file xray_reco_run189.h5 --runType=rtCalibration -b bGgPlot --ingrid --occupancy --config plotData.toml

which gives us the following plot:

With many more plots here: ./../Figs/statusAndProgress/xrayFingerRun/run189/

One very important plot: -> So the peak is at around 3 keV instead of about 8 keV, as the plot from Amptek in the section above pretends.

• [ ] Maybe at CAST they changed the target?

2.8.1. Window layout with correct window rotation

## Super dumb 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, chroma, unchained
proc hitsStrongback(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() =
  let black = color(0.0, 0.0, 0.0)
  let nPoints = 256
  var xs = linspace(-7.0, 7.0, nPoints)
  var dataX = newSeqOfCap[float](nPoints^2)
  var dataY = newSeqOfCap[float](nPoints^2)
  var inside = newSeqOfCap[bool](nPoints^2)
  for x in xs:
    for y in xs:
      if x*x + y*y < 7.0 * 7.0:
        when false:
          dataX.add x * cos(30.°.to(Radian)) + y * sin(30.°.to(Radian))
          dataY.add y * cos(30.°.to(Radian)) - x * sin(30.°.to(Radian))
          inside.add hitsStrongback(y)
        else:
          dataX.add x
          dataY.add y
          # rotate current y back, such that we can analyze in a "non rotated" coord. syst
          let yRot = y * cos(-30.°.to(Radian)) - x * sin(-30.°.to(Radian))
          inside.add hitsStrongback(yRot)
  let df = toDf(dataX, dataY, inside)
  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]") +
    xlim(-7, 7) + ylim(-7, 7) + 
    ggsave("/home/basti/org/Figs/statusAndProgress/detector/SiN_window_occlusion_rotated.png", width = 1150, height = 1000)
sample()

2.9.1. 2017/18 (Run 2)

Number of total events:                           3758960
Number of events without center:                  1557934 | 41.44587864728542%
Number of events only center:                     23820 | 0.633685913124907%
Number of events with center activity and outer:  984319 | 26.185939728009878%
Number of events any hit events:                  2542253 | 67.6318183752953%
Mean of event durations:                          2.144074329358038

Interestingly, the histogram of event durations looks as follows, fig. 6.

We can cut to the range between 0 and 2.2 s, fig. 7.

The peak at 0 is plain and simply a peak at exact 0 values (the previous figure only removed exact 0 values).

What does the energy distribution look like for these events? Fig. 8.

And the same split up per run (to make sure it's not one bad run), fig. 9.

Hmm. I suppose it's a bug in the firmware that the event duration is not correctly returned? Could happen if FADC triggers and for some reason 0 clock cycles are returned. This could be connected to the weird "hiccups" the readout sometimes does (when the FADC doesn't actually trigger for a full event). Maybe these are the events right after?

1. Noisy pixels
   

   In this run there are a few noisy pixels that need to be removed before background rates are calculated. These are listed in tab. 2.

   Table 2: Number of counts noisy pixels in 2017/18 dataset contribute to the number of background clusters remaining. The total number of noise clusters amounts to 1265 in this case (depends on the clustering algorithm potentially). These must be removed for a sane background level (and the area must be removed from from the size of active area in this dataset. NOTE: When using these numbers, make sure the x and y coordinates are not accidentally inverted.

   x	y	Count after logL
   64	109	7
   64	110	9
   65	108	30
   66	108	50
   67	108	33
   65	109	74
   66	109	262
   67	109	136
   68	109	29
   65	110	90
   66	110	280
   67	110	139
   65	111	24
   66	111	60
   67	111	34
   67	112	8

   \clearpage

2.9.2. End of 2018 (Run 3)

NOTE: <2022-12-20 Tue 15:30> In Run 3 we only used 2/30 as run settings! Hence a single peak in event duration.

And the same plots and numbers for 2018.

Number of total events:                           1837330
Number of events without center:                  741199 | 40.34109278137297%
Number of events only center:                     9462 | 0.514986420512374%
Number of events with center activity and outer:  470188 | 25.590830172043127%
Number of events any hit events:                  1211387 | 65.9319229534161%
Mean of event durations:                          2.1157526632342307

2.11.1. Information (raytracing, effective area etc) from CAST Nature paper

import ggplotnim
const path = "/home/basti/org/resources/llnl_cast_nature_jaime_data/2016_DEC_Final_CAST_XRT/EffectiveArea.txt"
const pathJCAP = "/home/basti/org/resources/llnl_xray_telescope_cast_effective_area.csv"

let dfJcap = readCsv(pathJCAP)
let df = readCsv(path, sep = ' ')
  .rename(f{"Energy[keV]" <- "E(keV)"},
          f{"EffectiveArea[cm²]" <- "Area(cm^2)"})
  .select("Energy[keV]", "EffectiveArea[cm²]")
let dfC = bind_rows([("JCAP", dfJcap), ("Nature", df)], "Type")
ggplot(dfC, aes("Energy[keV]", "EffectiveArea[cm²]", color = "Type")) +
  geom_line() +
  ggsave("/tmp/effective_area_jcap_vs_nature_llnl.pdf")

1. LLNL raytracing for axion image and CoolX X-ray finger
   

   The DTU thesis contains raytracing images (from page 78) for the X-ray finger run and for the axion image.

   1. X-ray finger
      

      The image (as a screenshot) from the X-ray finger:

      

      where we can see a few things:

      1. the caption mentions the source was 14.2 m away from the optic. This is nonsensical. The magnet is 9.26m long and even with the cryo housing etc. we won't get to much more than 10 m from the telescope. The X-ray finger was installed in the bore of the magnet!
      2. it mentions the source being 6 mm diameter (text mentions diameter explicitly). All we know about it is from the manufacturer that the size is given as 15 mm. But there is nothing about the actual size of the emission surface.
      3. the resulting raytraced image has a size of only slightly less than 3 mm in the short axis and maybe about 3 mm in the long axis.

      About 3: Our own X-ray finger is the following: file:///home/basti/phd/Figs/CAST_Alignment/xray_finger_centers_run_189.pdf (Note: it needs to be rotated of course) We can see that our real image is much larger! Along "x" it goes from about 5.5 to 10 mm or so! Quite a bit larger. And along y from less than 4 to maybe 10!

      Given that we have the raytracing data from Jaime, let's plot their data to see if it actually looks like that:

      import ggplotnim, sequtils, seqmath
      let df = readCsv("/home/basti/org/resources/llnl_cast_nature_jaime_data/2016_DEC_Final_CAST_XRT/3.00keV_2Dmap_CoolX.txt",
                       sep = ' ', skipLines = 2, colNames = @["x", "y", "z"])
        .mutate(f{"x" ~ `x` - mean(`x`)},
                f{"y" ~ `y` - mean(`y`)})
      var customInferno = inferno()
      customInferno.colors[0] = 0 # transparent
      ggplot(df, aes("x", "y", fill = "z")) +
        geom_raster() +
        scale_fill_gradient(customInferno) +
        xlab("x [mm]") + ylab("y [mm]") +
        ggtitle("LLNL raytracing of X-ray finger (Jaime)") +
        ggsave("~/org/Figs/statusAndProgress/rayTracing/raytracing_xray_finger_llnl_jaime.pdf")
      ggplot(df.filter(f{`x` >= -7.0 and `x` <= 7.0 and `y` >= -7.0 and `y` <= 7.0}),
             aes("x", "y", fill = "z")) +
        geom_raster() +
        scale_fill_gradient(customInferno) +
        xlab("x [mm]") + ylab("y [mm]") +
        xlim(-7.0, 7.0) + ylim(-7.0, 7.0) + 
        ggtitle("LLNL raytracing of X-ray finger zoomed (Jaime)") +
        ggsave("~/org/Figs/statusAndProgress/rayTracing/raytracing_xray_finger_llnl_jaime_gridpix_size.pdf")  
      

      This yields the following figure: and cropped to the range of the GridPix:

      This is MUCH bigger than the plot from the paper indicates. And the shape is also much more elongated! More in line with what we really see.

      Let's use our raytracer to produce the X-ray finger according to the specification of 14.2 m first and then a more reasonable estimate.

      Make sure to put the following into the config.toml file:

      [TestXraySource]
      useConfig = true # sets whether to read these values here. Can be overriden here or useng flag `--testXray`
      active = true     # whether the source is active (i.e. Sun or source?)
      sourceKind = "classical" # whether a "classical" source or the "sun" (Sun only for position *not* for energy)
      parallel = false
      energy = 3.0      # keV The energy of the X-ray source
      distance = 14200 # 9260.0 #106820.0 #926000 #14200 #9260.0 #2000.0 # mm  Distance of the X-ray source from the readout
      radius = 3.0 #21.5 #44.661 #8.29729 #46.609 #4.04043 #3.0 #4.04043 #21.5 # #21.5    # mm  Radius of the X-ray source
      offAxisUp = 0.0   # mm
      offAxisLeft = 0.0 # mm
      activity = 0.125  # GBq The activity in `GBq` of the source
      lengthCol = 0.0 #0.021 # mm  Length of a collimator in front of the source
      

      ./raytracer --ignoreDetWindow --ignoreGasAbs --suffix "_xrayFinger_14.2m_3mm"
      

      

      which more or less matches the size of our real data.

      Now the same with a source that is 10 m away:

      [TestXraySource]
      useConfig = true # sets whether to read these values here. Can be overriden here or useng flag `--testXray`
      active = true     # whether the source is active (i.e. Sun or source?)
      sourceKind = "classical" # whether a "classical" source or the "sun" (Sun only for position *not* for energy)
      parallel = false
      energy = 3.0      # keV The energy of the X-ray source
      distance = 10000 # 9260.0 #106820.0 #926000 #14200 #9260.0 #2000.0 # mm  Distance of the X-ray source from the readout
      radius = 3.0 #21.5 #44.661 #8.29729 #46.609 #4.04043 #3.0 #4.04043 #21.5 # #21.5    # mm  Radius of the X-ray source
      offAxisUp = 0.0   # mm
      offAxisLeft = 0.0 # mm
      activity = 0.125  # GBq The activity in `GBq` of the source
      lengthCol = 0.0 #0.021 # mm  Length of a collimator in front of the source
      

      ./raytracer --ignoreDetWindow --ignoreGasAbs --suffix "_xrayFinger_10m_3mm"
      

      

      which is quite a bit bigger than our real data. Maybe we allow some angles that we shouldn't, i.e. the X-ray finger has a collimator? Or our reflectivities are too good for too large angles?

      Without good knowledge of the real size of the X-ray finger emission this is hard to get right.

   2. Axion image
      

      The axion image as mentioned in the PhD thesis is the following:

      

      First of all let's note that the caption talks about emission of a 3 arcminute source. Let's check the apparent size of the sun and the typical emission, which is from the inner 30%:

      import unchained, math
      
      let Rsun = 696_342.km # SOHO mission 2003 & 2006
      # use the tangent to compute based on radius of sun:
      # tan α = Rsun / 1.AU
      echo "Apparent size of the sun = ", arctan(Rsun / 1.AU).Radian.to(ArcMinute)
      echo "Typical emission sun from inner 30% = ", arctan(Rsun * 0.3 / 1.AU).Radian.to(ArcMinute)
      let R3arc = (tan(3.ArcMinute.to(Radian)) * 1.AU).to(km)
      echo "Used radius for 3' = ", R3arc
      echo "As fraction of solar radius = ", R3arc / RSun
      

      So 3' correspond to about 18.7% of the radius. All in all that seems reasonable at least.

      Let's plot the axion image as we have it from Jaime's data:

      import ggplotnim, seqmath
      import std / [os, sequtils, strutils]
      proc readRT(p: string): DataFrame =
        result = readCsv(p, sep = ' ', skipLines = 4, colNames = @["x", "y", "z"])
        result["File"] = p  
      proc meanData(df: DataFrame): DataFrame =
        result = df.mutate(f{"x" ~ `x` - mean(col("x"))},
                           f{"y" ~ `y` - mean(col("y"))})
      
      proc plots(df: DataFrame, title, outfile: string) =
        var customInferno = inferno()
        customInferno.colors[0] = 0 # transparent
        ggplot(df, aes("x", "y", fill = "z")) +
          geom_raster() +
          scale_fill_gradient(customInferno) +
          xlab("x [mm]") + ylab("y [mm]") +
          ggtitle(title) +
          ggsave(outfile)
        ggplot(df.filter(f{`x` >= -7.0 and `x` <= 7.0 and `y` >= -7.0 and `y` <= 7.0}),
               aes("x", "y", fill = "z")) +
          geom_raster() +
          scale_fill_gradient(customInferno) +
          xlab("x [mm]") + ylab("y [mm]") +
          xlim(-7.0, 7.0) + ylim(-7.0, 7.0) + 
          ggtitle(title & " (zoomed)") +
          ggsave(outfile.replace(".pdf", "_gridpix_size.pdf"))
        
      block Single:    
        let df = readRT("/home/basti/org/resources/llnl_cast_nature_jaime_data/2016_DEC_Final_CAST_XRT/3.00keV_2Dmap.txt")
          .meanData()
        df.plots("LLNL raytracing of axion image @ 3 keV (Jaime)",
                 "~/org/Figs/statusAndProgress/rayTracing/raytracing_axion_image_llnl_jaime_3keV.pdf")
      
      block All:
        var dfs = newSeq[DataFrame]()
        for f in walkFiles("/home/basti/org/resources/llnl_cast_nature_jaime_data/2016_DEC_Final_CAST_XRT/*2Dmap.txt"):
          echo "Reading: ", f
          dfs.add readRT(f)
        echo "Summarize"
        var df = dfs.assignStack()
        df = df.group_by(@["x", "y"])
          .summarize(f{float: "z" << sum(`z`)},
                     f{float: "zMean" << mean(`z`)})
        df.writeCsv("/tmp/llnl_raytracing_jaime_all_energies_raw_sum.csv")
        df = df.meanData()
        df.writeCsv("/tmp/llnl_raytracing_jaime_all_energies.csv")  
        plots(df,
              "LLNL raytracing of axion image (sum all energies) (Jaime)",
              "~/org/Figs/statusAndProgress/rayTracing/raytracing_axion_image_llnl_jaime_all_energies.pdf")
      

      The 3 keV data for the axion image: and cropped again:

      And the sum of all energies: and cropped again:

      Both clearly show the symmetric shape that is so weird but also - again - does NOT reproduce the raytracing seen in the screenshot above! That one clearly has a very stark tiny center with the majority of the flux, which is gone and replaced by a much wider region of significant flux!

      Both are in strong contrast to our own axion image. Let's compute that using the Primakoff only (make sure to disable the X-ray test source in the config file!):

      ./raytracer --ignoreDetWindow --ignoreGasAbs --suffix "_axionImagePrimakoff_focal_point"
      

      

      and for a more realistic image at the expected conversion point:

      [DetectorInstallation]
      useConfig = true # sets whether to read these values here. Can be overriden here or using flag `--detectorInstall`
      # Note: 1500mm is LLNL focal length. That corresponds to center of the chamber!
      distanceDetectorXRT = 1487.93  # mm
      distanceWindowFocalPlane = 0.0 # mm
      lateralShift = 0.0             # mm lateral ofset of the detector in repect to the beamline
      transversalShift = 0.0         # mm transversal ofset of the detector in repect to the beamline #0.0.mm #
      

      ./raytracer --ignoreDetWindow --ignoreGasAbs --suffix "_axionImagePrimakoff_conversion_point"
      

      which yields:

      

      which is not that far off in size of the LLNL raytraced image. The shape is just quite different!

2. Reply to Igor about LLNL telescope raytracing
   

   Igor wrote me the following mail:

   Hi Sebastian, Now that we are checking with Cristina the shape of the signal after the LLNL telescope for the SRMM analysis, I got two questions on your analysis:

   1. The signal spot shape that you present is different from the one we have for the Nature physics paper. Do you understand why? There was a change in the Ingrid setup wrt the SRMM setup that explains it, maybe?
   2. Do you have a spot calibration data that allows to crosscheck the position (and rotation) of the signal spot in the Ingrid chip coordinates?

   Best, Igor

   as a reply to my "Limit method for 7-GridPix @ CAST" mail on <2023-08-09 Wed>.

   I ended up writing a lengthy reply.

   The reply is also found here: ./../Mails/igorReplyLLNL/igor_reply_llnl_axion_image.html

   1. My reply
      

      Hey,

      sorry for the late reply. I didn't want to reply with one sentence for each question. While looking into the questions in more details more things came up.

      One thing - embarrassingly - is that I completely forgot to apply the rotation of my detector in the limit calculation (in our case the detector is rotated by 90° compared to the "data" x-y plane). Added to that is the slight rotation of the LLNL axis, which I also need to include (here I simply forgot that we never added it to the raytracer. Given that the spacer is not visible in the axion image, it didn't occur to me).

      Let's start with your second question

      Do you have a spot calibration data that allows to crosscheck the position (and rotation) of the signal spot in the Ingrid chip coordinates?

      Yes, we have two X-ray finger runs. Unfortunately, one of them is not useful, as it was taken in July 2017 after our detector had to be removed again to make space for a short KWISP data taking. We have a second one from April 2018, which is partially useful. However, the detector was again dismounted between April and October 2018 and we don't have an X-ray finger run for the last data taking between Oct 2018 to Dec 2018.

      Fig. 14 shows the X-ray finger run shows the latter X-ray finger run. The two parallel lines with few clusters are two of the window strongbacks. The other line is the graphite spacer of the telescope. The center positions of the clusters are at

      • (x, y) = (7.43, 6.59)

      (the chip center is at (7, 7). This is what makes up the basis of our position systematic uncertainty of 5%. The 5% correspond to 0.05*7 mm = 0.35 mm.

      

      I decided not to move the actual center of the solar axion image because the X-ray finger data is hard to interpret for three different reasons:

      1. The entire CAST setup is "modified" in between normal data takings and installation of the X-ray finger. Who knows the effect warming up the magnet etc. is on the spot position?
      2. determining the actual center position of the axion spot based on the X-ray finger cluster centers is problematic due to the fact that the LLNL telescope is only a portion of a full telescope. With the resulting shape of the X-ray finger signal, combined with the missing data due to the window strongback and graphite spacer and relatively low statistics in the first place, makes trusting the numbers problematic.
      3. as I said before, we don't even have an X-ray finger run for the last part of the data taking. While we have the geometer measurements from the targets, I don't have the patience to learn about the coordinate system they use and attempt to reconstruct the possible movement based on those measured coordinates.

      Given that we take into account the possible movement in the systematics, I believe this is acceptable.

      The signal spot shape that you present is different from the one we have for the Nature physics paper. Do you understand why? There was a change in the Ingrid setup wrt the SRMM setup that explains it, maybe?

      Here we now come to the actual part that is frustrating for me, too. Unfortunately, due to the "black box" nature of the LLNL telescope, Johanna and me never fully understood this. We don't understand how the raytracing calculations done by Michael Pivovaroff can ever produce a symmetric image given that the LLNL telescope is a) not a perfect Wolter design, but has cone shaped mirrors, b) is only a small portion of a full telescope and c) the incoming X-rays are not perfectly parallel. Intuitively I don't expect to have a symmetric image there. And our raytracing result does not produce anything like that.

      A couple of years ago Johanna tried to find out more information about the LLNL raytracing results, but back then when Julia and Jaime were still at LLNL, the answer was effectively a "it's a secret, we can't provide more information".

      As such all I can do is try to reproduce the results as well as possible. If they don't agree all I can do is provide explanations about what we compute and give other people access to my data, code and results. Then at least we can all hopefully figure out if there's something wrong with our approach.

      Fig. 15 is the raytracing result as it is presented on page 78 of the PhD thesis of A. Jakobsen. It mentions that the Sun is considered as a 3' source, implying the inner ~18% of the Sun are contributing to axion emission.

      

      If I compute this with our own raytracer for the focal spot, I get the plot shown in fig. \ref{fig:axion_image_primakoff_focal_spot}. Fig. \ref{fig:axion_image_primakoff_median_conv} then corresponds to the point that sees the median of all conversions in the gas based on X-ray absorption in the gas. This is now for the case of a pure Primakoff emission and not for dominant axion-electron coupling, as I showed in my presentation (this changes the dominant contributions by radius slightly, see fig. \ref{fig:radial_production_primakoff} { Primakoff } and fig. \ref{fig:radial_production_axion_electron} { axion-electron }). They look very similar, but there are slight changes between the two axion images.

      This is one of the big reasons I want to have my own raytracing simulation. Different emission models result in different axion images!

      \begin{figure}[htbp] \centering \begin{subfigure}{0.5\linewidth} \centering \includegraphics[width=0.95\textwidth]{/home/basti/org/Figs/statusAndProgress/rayTracing/raytracing_axion_image_primakoff_focal_point.pdf} \caption{Focal spot} \label{fig:axion_image_primakoff_focal_spot} \end{subfigure}% \begin{subfigure}{0.5\linewidth} \centering \includegraphics[width=0.95\textwidth]{/home/basti/org/Figs/statusAndProgress/rayTracing/raytracing_axion_image_primakoff_conversion_point.pdf} \caption{Median conversion point} \label{fig:axion_image_primakoff_median_conv} \end{subfigure} \label{fig:axion_image} \caption{\subref{fig:axion_image_primakoff_focal_spot} Axion image for Primakoff emission from the Sun, computed for the exact LLNL focal spot. (Ignore the title) \subref{fig:axion_image_primakoff_median_conv} Axion image for the median conversion point of the X-rays actually entering the detector. } \end{figure} \begin{figure}[htbp] \centering \begin{subfigure}{0.5\linewidth} \centering \includegraphics[width=0.95\textwidth]{~/org/Figs/statusAndProgress/axionProduction/sampled_radii_primakoff.pdf} \caption{Primakoff radii} \label{fig:radial_production_primakoff} \end{subfigure}% \begin{subfigure}{0.5\linewidth} \centering \includegraphics[width=0.95\textwidth]{~/org/Figs/statusAndProgress/axionProduction/sampled_radii_axion_electron.pdf} \caption{Axion-electron radii} \label{fig:radial_production_axion_electron} \end{subfigure} \label{fig:radial_production} \caption{\subref{fig:radial_production_primakoff} Radial production in the Sun for Primakoff emission. \subref{fig:radial_production_axion_electron} Radial production for axion-electron emission. } \end{figure}

      Note that this currently does not yet take into account the slight rotation of the telescope. I first need to extract the rotation angle from the X-ray finger run.

      Fig. 16 is the sum of all energies of the raytracing results that Jaime finally sent to Cristina a couple of weeks ago. In this case cropped to the size of our detector, placed at the center. These should be - as far as I understand - the ones that the contours used in the Nature paper are based on. However, these clearly do not match the results shown in the PhD thesis of Jakobsen. The extremely small focus area in black is gone and replaced by a much more diffuse area. But again, it is very symmetric, which I don't understand.

      

      And while I was looking into this I also thought I should try to (attempt to) reproduce the X-ray finger raytracing result. Here came another confusion, because the raytracing results for that shown in the PhD thesis, fig. 17, mention that the X-ray finger was placed \SI{14.2}{m} away from the optic with a diameter of \SI{6}{mm}. That seems very wrong, given that the magnet bore is only \SI{9.26}{m} long. In total the entire magnet is - what - maybe \SI{10}{m}? At most it's maybe \SI{11}{m} to the telescope when the X-ray finger is installed in the bore? Unfortunately, the website about the X-ray finger from Amptek is not very helpful either:

      https://www.amptek.com/internal-products/obsolete-products/cool-x-pyroelectric-x-ray-generator

      as the only thing it says about the size is:

      Miniature size: 0.6 in dia x 0.4 in (15 mm dia x 10 mm)

      

      Nothing about the actual size of the area that emits X-rays. Neither do I know anything about a possible collimator used.

      Furthermore, the spot size seen here is only about \(\sim 2.5·\SI{3}{mm²}\) or so. Comparing it to the spot size seen with our detector it's closer to \(\sim 5·\SI{5}{mm²}\) or even a bit larger!

      So I decided to run a raytracing following these numbers, i.e. \(\SI{14.2}{m}\) and a \(\SI{3}{mm}\) radius disk shaped source without a collimator. That yields fig. 18. As we can see the size is more in line with our actually measured data.

      

      Again, I looked at the raytracing results that Jaime sent to Cristina, which includes a file with suffix "CoolX". That plot is shown in fig. 19. As we can see, it is also much larger suddenly than shown in the PhD thesis (more than \(4 · \SI{4}{mm²}\)), slightly smaller than ours.

      Note that the Nature paper mentions the source is about \(\SI{12}{m}\) away. I was never around when the X-ray finger was installed, nor do I have any good data about the real magnet size or lengths of the pipes between magnet and telescope.

      

      So, uhh, yeah. This is all very confusing. No matter where one looks regarding this telescope, one is bound to find contradictions or just confusing statements… :)

2.11.2. Information from NuSTAR PhD thesis

2.11.3. X-ray raytracers

2.11.4. DTU FTP server [/]

2.11.5. Michael Pivovaroff talk about Axions, CAST, IAXO

3.1.1. Solar axion-electron flux

import ggplotnim
const path = "~/org/resources/differential_flux_sun_earth_distance/solar_axion_flux_differential_g_ae_1e-13_g_ag_1e-12_g_aN_1e-15_0.989AU.csv"
let df = readCsv(path)
  .filter(f{`type` notin ["LP Flux", "TP Flux", "57Fe Flux"]})
echo df
ggplot(df, aes("Energy", "diffFlux", color = "type")) +
  geom_line() +
  xlab(r"Energy [$\si{keV}$]", margin = 1.5) + ylab(r"Flux [$\si{keV^{-1}.cm^{-2}.s^{-1}}$]", margin = 2.75) +
  ggtitle(r"Differential solar axion flux for $g_{ae} = \num{1e-13}, g_{aγ} = \SI{1e-12}{GeV^{-1}}$") +
  xlim(0, 10) +
  margin(top = 1.5, left = 3.25) +
  theme_transparent() + 
  ggsave("~/org/Figs/statusAndProgress/differential_flux_sun_earth_distance/differential_solar_axion_fluxg_ae_1e-13_g_ag_1e-12_g_aN_1e-15_0.989AU.pdf", width = 800, height = 480,
         useTeX = true, standalone = true)        

3.1.2. Radial production

solarModelFile  = "solar_model_dataframe.csv"

./raytracer --ignoreDetWindow --ignoreGasAbs --suffix "_axion_electron" --sanity

solarModelFile  = "solar_model_dataframe_fluxKind_fkAxionPhoton_0.989AU.csv" #solar_model_dataframe.csv"

./raytracer --ignoreDetWindow --ignoreGasAbs --suffix "_primakoff" --sanity

3.2.1. Deriving the missing constants in the conversion probability

3.3.1. Average distance X-rays travel in Argon at CAST conditions [/]

1. Analytical approach
   

   import xrayAttenuation, ggplotnim, unchained
   
   # 1. read the file containing efficiencies
   var effDf = readCsv("/home/basti/org/resources/combined_detector_efficiencies.csv")
     .mutate(f{"NoGasEff" ~ idx("300nm SiN") * idx("20nm Al") * `LLNL`})
   
   # 2. compute the absorption length for Argon
   let ar = Argon.init()
   let ρ_Ar = density(1050.mbar.to(Pascal), 293.K, ar.molarMass)
   effDf = effDf
     .filter(f{idx("Energy [keV]") > 0.05})
     .mutate(f{float: "l_abs" ~ absorptionLength(ar, ρ_Ar, idx("Energy [keV]").keV).float})
   # compute the weighted mean of the effective flux behind the window with the
   # absorption length, i.e.
   # `<x> = Σ_i (ω_i x_i) / Σ_i ω_i`
   let weightedMean = (effDf["NoGasEff", float] *. effDf["l_abs", float]).sum() /
     effDf["NoGasEff", float].sum()
   echo "Weighted mean of distance: ", weightedMean.Meter.to(cm)
   # for reference the effective flux:
   ggplot(effDf, aes("Energy [keV]", "NoGasEff")) +
     geom_line() +
     ggsave("/tmp/combined_efficiency_no_gas.pdf")
   ggplot(effDf, aes("Energy [keV]", "l_abs")) +
     geom_line() +
     ggsave("/tmp/absorption_length_argon_cast.pdf")
   

   This means the "effective" position of the axion image should be 0.0122 m or 1.22 cm in the detector. This is (fortunately) relatively close to the 1.5 cm (center of the detector) that we used so far.

   • [X]

     Is the above even correct? The absorption length describes the distance at which only 1/e particles are left. That means at that distance (1 - 1/e) have disappeared. To get a number don't we need to do a monte carlo (or some kind of integral) of the average? -> Well, the mean of an exponential distribution is 1/λ (if defined as \(\exp(-λx)\)!), from that point of view I think the above is perfectly adequate! Note however that the median of the distribution is \(\frac{\ln 2}{λ}\)! When looking at the distribution of our transverse RMS values for example the peak corresponds to something that is closer to the median (but is not exactly the median either; the peak is the 'mode' of the distribution). Arguably more interesting is the cutoff we see in the data as that corresponds to the largest possible diffusion (but again that is being folded with the statistics of getting a larger RMS! :/ )

     UPDATE: <2023-06-14 Wed 17:24> See the section below for the numerical approach. As it turns out the above unfortunately is not correct for 3 important reasons (2 of which we were aware of):

     1. It does not include the axion spectrum, it changes the location of the mean slightly.
     2. It implicitly assumes all X-rays of all energies will be detected. This implies an infinitely long detector and not our detector limited by a length of 3 cm! This skews the actual mean to lower values, because the mean of those that are detected are at smaller values.
     3. Point 2 implies not only that some X-rays won't be detected, but effectively it gives a higher weight to energies that are absorbed with certainty compared to those that sometimes are not absorbed! This further reduces the mean. This can be interpreted as reducing the input flux by the percentage of the absorption probability for each energy. In this sense the above needs to be multiplied by the absorption probability to be more correct! Yet this still does not make it completely right, as that just assumes the fraction of photons of a given energy are reduced, but not that all detected ones have lengths consistent with a 3cm long volume!
     4. (minor) does not include isobutane.

   A (shortened and) improved version of the above (but still not quite correct!):

   import xrayAttenuation, ggplotnim, unchained
   # 1. read the file containing efficiencies
   var effDf = readCsv("/home/basti/org/resources/combined_detector_efficiencies.csv")
     .mutate(f{"NoGasEff" ~ idx("300nm SiN") * idx("20nm Al") * `LLNL` * idx("30mm Ar Abs.")})
   # 2. compute the absorption length for Argon
   let ar = Argon.init()
   let ρ_Ar = density(1050.mbar.to(Pascal), 293.K, ar.molarMass)
   effDf = effDf.filter(f{idx("Energy [keV]") > 0.05})
     .mutate(f{float: "l_abs" ~ absorptionLength(ar, ρ_Ar, idx("Energy [keV]").keV).float})
   let weightedMean = (effDf["NoGasEff", float] *. effDf["l_abs", float]).sum() /
     effDf["NoGasEff", float].sum()
   echo "Weighted mean of distance: ", weightedMean.Meter.to(cm)
   

   We could further multiply in the axion flux of course, but as this cannot be fully correct in this way, we'll do it numerically. We would have to calculate the real mean of the exponential distribution for each energy based on the truncated exponential distribution. Effectively we have a bonded exponential between 0 and 3 cm, whose mean is of course going to differ from the parameter \(λ\).

2. Numerical approach
   

   Let's write a version of the above code that computes the result by sampling from the exponential distribution for the conversion point.

   What we need:

   • our sampling logic
   • sampling from exponential distribution depending on energy
   • the axion flux

   Let's start by importing the modules we need:

   import helpers / sampling_helper # sampling distributions
   import unchained                 # sane units
   import ggplotnim                 # see something!
   import xrayAttenuation           # window efficiencies
   import math, sequtils
   

   where the sampling_helpers is a small module to sample from a procedure or a sequence.

   In addition let's define some helpers:

   from os import `/`
   const ResourcePath = "/home/basti/org/resources"
   const OutputPath = "/home/basti/org/Figs/statusAndProgress/axion_conversion_point_sampling/"
   

   Now let's read the LLNL telescope efficiency as well as the axion flux model. Note that we may wish to calculate the absorption points not only for a specific axion flux model, but potentially any other kind of signal. We'll build in functionality to disable different contributions.

   let dfAx = readCsv(ResourcePath / "solar_axion_flux_differential_g_ae_1e-13_g_ag_1e-12_g_aN_1e-15.csv")
     .filter(f{`type` == "Total flux"})
   let dfLLNL = readCsv(ResourcePath / "llnl_xray_telescope_cast_effective_area_parallel_light_DTU_thesis.csv")
     .mutate(f{"Efficiency" ~ idx("EffectiveArea[cm²]") / (PI * 2.15 * 2.15)})
   

   Note: to get the differential axion flux use readOpacityFile from https://github.com/jovoy/AxionElectronLimit. It generates the CSV file.

   Next up we need to define the material properties of the detector window in order to compute its transmission.

   let Si₃N₄ = compound((Si, 3), (N, 4)) # actual window
   const ρSiN = 3.44.g•cm⁻³
   const lSiN = 300.nm                  # window thickness
   let Al = Aluminium.init()            # aluminium coating
   const ρAl = 2.7.g•cm⁻³
   const lAl = 20.nm                    # coating thickness
   

   With these numbers we can compute the transmission at an arbitrary energy. In order to compute the correct inputs for the calculation we now have everything. We wish to compute the following, the intensity \(I(E)\) is the flux that enters the detector

   \[ I(E) = f(E) · ε_{\text{LLNL}} · ε_{\ce{Si3.N4}} · ε_{\ce{Al}} \]

   where \(f(E)\) is the solar axion flux and the \(ε_i\) are the efficiencies associated with the telescope and transmission of the window. The idea is to sample from this intensity distribution to get a realistic set of X-rays as they would be experienced in the experiment. One technical aspect still to be done is an interpolation of the axion flux and LLNL telescope efficiency to evaluate the data at an arbitrary energy as to define a function that yields \(I(E)\).

   Important note: We fully neglect here the conversion probability and area of the magnet bore. These (as well as a potential time component) are purely constants and do not affect the shape of the distribution \(I(E)\). We want to sample from it to get the correct weighting of the different energies, but do not care about absolute numbers. So differential fluxes are fine.

   The idea is to define the interpolators and then create a procedure that captures the previously defined properties and interpolators.

   from numericalnim import newLinear1D, eval
   let axInterp = newLinear1D(dfAx["Energy", float].toSeq1D, dfAx["diffFlux", float].toSeq1D)
   let llnlInterp = newLinear1D(dfLLNL["Energy[keV]", float].toSeq1D, dfLLNL["Efficiency", float].toSeq1D)
   

   With the interpolators defined let's write the implementation for \(I(E)\):

   proc I(E: keV): float =
     ## Compute the intensity of the axion flux after telescope & window eff.
     ##
     ## Axion flux and LLNL efficiency can be disabled by compiling with
     ## `-d:noAxionFlux` and `-d:noLLNL`, respectively.
     result = transmission(Si₃N₄, ρSiN, lSiN, E) * transmission(Al, ρAl, lAl, E)
     when not defined(noAxionFlux):
       result *= axInterp.eval(E.float)
     when not defined(noLLNL):
       result *= llnlInterp.eval(E.float)
     
   

   Let's test it and see what we get for e.g. \(\SI{1}{keV}\):

   echo I(1.keV)
   

   yields \(1.249e20\). Not the most insightful, but it seems to work. Let's plot it:

   let energies = linspace(0.01, 10.0, 1000).mapIt(it.keV)
   let Is = energies.mapIt(I(it))
   block PlotI:
     let df = toDf({ "E [keV]" : energies.mapIt(it.float),
                     "I" : Is })
     ggplot(df, aes("E [keV]", "I")) +
       geom_line() +
       ggtitle("Intensity entering the detector gas") + 
       ggsave(OutputPath / "intensity_axion_conversion_point_simulation.pdf")
   

   shown in fig. 20. It looks exactly as we would expect.

   

   Now we define the sampler for the intensity distribution \(I(E)\), which returns an energy weighted by \(I(E)\):

   let Isampler = sampler(
     (proc(x: float): float = I(x.keV)), # wrap `I(E)` to take `float`
     0.01, 10.0, num = 1000 # use 1000 points for EDF & sample in 0.01 to 10 keV
   )
   

   and define a random number generator:

   import random
   var rnd = initRand(0x42)
   

   First we will sample 100,000 energies from the distribution to see if we recover the intensity plot from before.

   block ISampled:
     const nmc = 100_000
     let df = toDf( {"E [keV]" : toSeq(0 ..< nmc).mapIt(rnd.sample(Isampler)) })
     ggplot(df, aes("E [keV]")) +
       geom_histogram(bins = 200, hdKind = hdOutline) +
       ggtitle("Energies sampled from I(E)") +
       ggsave(OutputPath / "energies_intensity_sampled.pdf")
   

   This yields fig. 21, which clearly shows the sampling works as intended.

   

   The final piece now is to use the same sampling logic to generate energies according to \(I(E)\), which correspond to X-rays of said energy entering the detector. For each of these energies then sample from the Beer-Lambert law

   \[ I(z) = I_0 \exp\left[ - \frac{z}{l_{\text{abs}} } \right] \] where \(I_0\) is some initial intensity and \(l_\text{abs}\) the absorption length. The absorption length is computed from the gas mixture properties of the gas used at CAST, namely Argon/Isobutane 97.7/2.3 at \(\SI{1050}{mbar}\). It is the inverse of the attenuation coefficient \(μ_M\)

   \[ l_{\text{abs}} = \frac{1}{μ_M} \]

   where the attenuation coefficient is computed via

   \[ μ_m = \frac{N_A}{M * σ_A} \]

   with \(N_A\) Avogadro's constant, \(M\) the molar mass of the compound and \(σ_A\) the atomic absorption cross section. The latter again is defined by

   \[ σ_A = 2 r_e λ f₂ \]

   with \(r_e\) the classical electron radius, \(λ\) the wavelength of the X-ray and \(f₂\) the second scattering factor. Scattering factors are tabulated for different elements, for example by NIST and Henke. For a further discussion of this see the README and implementation of xrayAttenuation.

   We will now go ahead and define the CAST gas mixture:

   proc initCASTGasMixture(): GasMixture =
     ## Returns the absorption length for the given energy in keV for CAST
     ## gas conditions:
     ## - Argon / Isobutane 97.7 / 2.3 %
     ## - 20°C ( for this difference in temperature barely matters)
     let arC = compound((Ar, 1)) # need Argon gas as a Compound
     let isobutane = compound((C, 4), (H, 10))
     # define the gas mixture
     result = initGasMixture(293.K, 1050.mbar, [(arC, 0.977), (isobutane, 0.023)])
   let gm = initCASTGasMixture()  
   

   To sample from the Beer-Lambert law with a given absorption length we also define a helper that returns a sampler for the target energy using the definition of a normalized exponential distribution

   \[ f_e(x, λ) = \frac{1}{λ} \exp \left[ -\frac{x}{λ} \right] \]

   The sampling of the conversion point is the crucial aspect of this. Naively we might want to sample between the detector volume from 0 to \(\SI{3}{cm}\). However, this skews our result. Our calculation depends on the energy distribution of the incoming X-rays. If the absorption length is long enough the probability of reaching the readout plane and thus not being detected is significant. Restricting the sampler to \(\SI{3}{cm}\) would pretend that independent of absorption length we would always convert within the volume, giving too large a weight to the energies that should sometimes not be detected!

   Let's define the sampler now. It takes the gas mixture and the target energy. A constant SampleTo is defined to adjust the position to which we sample at compile time (to play around with different numbers).

   proc generateSampler(gm: GasMixture, targetEnergy: keV): Sampler =
     ## Generate the exponential distribution to sample from based on the
     ## given absorption length
     # `xrayAttenuation` `absorptionLength` returns number in meter!
     let λ = absorptionLength(gm, targetEnergy).to(cm)
     let fnSample = (proc(x: float): float =
                       result = expFn(x, λ.float) # expFn = 1/λ · exp(-x/λ)
     )
     const SampleTo {.intdefine.} = 20 ## `SampleTo` can be set via `-d:SampleTo=<int>`
     let num = (SampleTo.float / 3.0 * 1000).round.int # number of points to sample at
     result = sampler(fnSample, 0.0, SampleTo, num = num)
   

   Note that this is inefficient, because we generate a new sampler from which we only sample a single point, namely the conversion point of that X-ray. If one intended to perform a more complex calculation or wanted to sample orders of magnitude more X-rays, one should either restructure the code (i.e. sample from known energies and then reorder based on the weight defined by \(I(E)\) or cache the samplers and pre-bin the energies.

   For reference let's compute the absorption length as a function of energy for the CAST gas mixture:

   block GasAbs:
     let df = toDf({ "E [keV]" : linspace(0.03, 10.0, 1000),
                     "l_abs [cm]" : linspace(0.03, 10.0, 1000).mapIt(absorptionLength(gm, it.keV).m.to(cm).float) })
     ggplot(df, aes("E [keV]", "l_abs [cm]")) +
       geom_line() +
       ggtitle("Absorption length of X-rays in CAST gas mixture: " & $gm) +
       margin(top = 1.5) + 
       ggsave(OutputPath / "cast_gas_absorption_length.pdf")
   

   which yields fig. 22

   

   So, finally: let's write the MC sampling!

   const nmc = 500_000 # start with 100k samples
   var Es = newSeqOfCap[keV](nmc)
   var zs = newSeqOfCap[cm](nmc)
   while zs.len < nmc:
     # 1. sample an energy according to `I(E)`
     let E = rnd.sample(Isampler).keV
     # 2. get the sampler for this energy
     let distSampler = generateSampler(gm, E) 
     # 3. sample from it
     var z = Inf.cm
     when defined(Equiv3cmSampling): ## To get the same result as directly sampling
                                     ## only up to 3 cm use the following code
       while z > 3.0.cm:
         z = rnd.sample(distSampler).cm 
     elif defined(UnboundedVolume): ## This branch pretends the detection volume
                                    ## is unbounded if we sample within 20cm
       z = rnd.sample(distSampler).cm 
     else: ## This branch is the physically correct one. If an X-ray reaches the
           ## readout plane it is _not_ recorded, but it was still part of the
           ## incoming flux!
       z = rnd.sample(distSampler).cm
       if z > 3.0.cm: continue # just drop this X-ray
     zs.add z
     Es.add E
   

   Great, now we have sampled the conversion points according to the correct intensity. We can now ask for statistics or create different plots (e.g. conversion point by energies etc.).

   import stats, seqmath # mean, variance and percentile
   let zsF = zs.mapIt(it.float) # for math
   echo "Mean conversion position = ", zsF.mean().cm
   echo "Median conversion position = ", zsF.percentile(50).cm
   echo "Variance of conversion position = ", zsF.variance().cm
   

   This prints the following:

   Mean conversion position = 0.556813 cm
   Median conversion position = 0.292802 cm
   Variance of conversion position = 0.424726 cm
   

   As we can see (unfortunately) our initial assumption of a mean distance of \(\SI{1.22}{cm}\) are quite of the mark. The more realistic number is only \(\SI{0.56}{cm}\). And if we were to use the median it's only \(\SI{0.29}{cm}\).

   Let's plot the conversion points of all sampled (and recorded!) X-rays as well as what their distribution against energy looks like.

   let dfZ = toDf({ "E [keV]" : Es.mapIt(it.float),
                    "z [cm]"  : zs.mapIt(it.float) })
   ggplot(dfZ, aes("z [cm]")) +
     geom_histogram(bins = 200, hdKind = hdOutline) +
     ggtitle("Conversion points of all sampled X-rays according to I(E)") +
     ggsave(OutputPath / "sampled_axion_conversion_points.pdf")
   ggplot(dfZ, aes("E [keV]", "z [cm]")) +
     geom_point(size = 1.0, alpha = 0.2) + 
     ggtitle("Conversion points of all sampled X-rays according to I(E) against their energy") +
     ggsave(OutputPath / "sampled_axion_conversion_points_vs_energy.png",
           width = 1200, height = 800)  
   

   The former is shown in fig. 23. The overlapping exponential distribution is obvious, as one would expect. The same data is shown in fig. 24, but in this case not as a histogram, but by their energy as a scatter plot. We can clearly see the impact of the absorption length on the conversion points for each energy!

   

   

3. Compiling and running the code
   

   The code above is written in literate programming style. To compile and run it we use ntangle to extract it from the Org file:

   ntangle <file>
   

   which generates ./../../../../tmp/sample_axion_xrays_conversion_points.nim.

   Compiling and running it can be done via:

   nim r -d:danger /tmp/sample_axion_xrays_conversion_points.nim
   

   which compiles and runs it as an optimized build.

   We have the following compilation flags to compute different cases:

   • -d:noLLNL: do not include the LLNL efficiency into the input intensity
   • -d:noAxionFlux: do not include the axion flux into the input intensity
   • -d:SampleTo=<int>: change to where we sample the position (only to 3cm for example)
   • -d:UnboundedVolume: if used together with the default SampleTo (or any large value) will effectively compute the case of an unbounded detection volume (i.e. every X-ray recorded with 100% certainty).
   • -d:Equiv3cmSampling: Running this with the default SampleTo (or any large value) will effectively change the sampling to a maximum \SI{3}{cm} sampling. This can be used as a good crossheck to verify the sampling behavior is independent of the sampling range.

   Configurations of note:

   nim r -d:danger -d:noAxionFlux /tmp/sample_axion_xrays_conversion_points.nim
   

   \(⇒\) realistic case for a flat input spectrum Yields:

   Mean conversion position = 0.712102 cm
   Median conversion position = 0.445233 cm
   Variance of conversion position = 0.528094 cm
   

   nim r -d:danger -d:noAxionFlux -d:UnboundedVolume /tmp/sample_axion_xrays_conversion_points.nim
   

   \(⇒\) the closest analogue to the analytical calculation from section 3.3.1.1 (outside of including isobutane here) Yields:

   Mean conversion position = 1.25789 cm
   Median conversion position = 0.560379 cm
   Variance of conversion position = 3.63818 cm
   

   nim r -d:danger /tmp/sample_axion_xrays_conversion_points.nim
   

   \(⇒\) the case we most care about and of which the numbers are mentioned in the text above.

4. Absorption edge in data
   

   Question:

   • [X] Can we see the absorption edge of Argon in our data? E.g. in the transverse RMS of the CDL data? In theory we should see a huge jump in the transverse nature (and cluster size) of the clusters above and below that point. MAYBE this could also relate to the strong cutoff we see in our background rate at 3 keV due to some effect of the efficiency of our cuts changing significantly there?

   If my "theory" is correct it would mean that the transverse RMS should be significantly different if I cut to the energy for e.g. the photo peak and escape peak?

   Update <2023-03-29 Wed 17:10>: As explained in multiple places since the above two TODOs were written. It's not as straightforward, because the exponential distribution still implies that a large fraction of events convert close to the cathode. The result is a smoothed out distribution of the RMS data, making the difference between escape and photo peak for example not as extreme as one might imagine. See the simulations below and the related FADC rise time simulations for more insight.

3.3.2. Simulating longitudinal and transverse cluster sizes using MC

import std / [random, sequtils, algorithm]
import seqmath, ggplotnim
template toEDF(data: seq[float], isCumSum = false): untyped =
  ## Computes the EDF of binned data
  var dataCdf = data
  if not isCumSum:
    seqmath.cumsum(dataCdf)
  let integral = dataCdf[^1]
  let baseline = min(data) # 0.0
  dataCdf.mapIt((it - baseline) / (integral - baseline))

proc sample(cdf: seq[float], ys: seq[float]): float =
  let point = rand(1.0)
  let idx = cdf.lowerBound(point)
  if idx < cdf.len:
    result = ys[idx]
  else:
    result = Inf

proc expFn(x: float, λ: float): float =
  result = 1.0 / λ * exp(- x / λ)

const Upper = 3.0
const λ = 2.0

let xs = linspace(0.0, Upper, 1000)
let ys = xs.mapIt(expFn(it, λ))

# now sample 100,000 points
let cdf = ys.toEdf()
let ySampled = toSeq(0 ..< 1_000_000).mapIt(sample(cdf, xs))
let dfS = toDf(ySampled)

ggplot(toDf(xs, cdf), aes("xs", "cdf")) +
  geom_line() +
  ggsave("/t/test_cdf.pdf")

echo dfS

# rescale according to normalization of the range we use
# normalize by y = y / (∫_Lower^Upper f(x) dx) =  # Lower = 0, Upper = 3.0 (`Upper`)
#              y = y / (∫_0^Upper 1/λ exp(-x/λ) dx = y / [ ( -exp(-x/λ) )|^Upper_0 ] 
#              y = y / [ (-exp(-Upper/λ) - (-exp(-Lower/λ) ) ]
#              y = y / [ (-exp(-3.0/λ)) + 1 ]   ^--- 1 = exp(0)
let df = toDf(xs, ys)
  .mutate(f{"ys" ~ `ys` / (-exp(-Upper / λ) + 1.0)})

ggplot(df, aes("xs")) +
  geom_line(aes = aes(y = "ys")) +
  geom_histogram(data = dfS, aes = aes(x = "ySampled"), bins = 100, density = true, alpha = 0.5, hdKind = hdOutline, fillColor = "red") +
  ggsave("/t/test_sample.pdf")

import std / [random, sequtils, algorithm]
import seqmath, ggplotnim
template toEDF*(data: seq[float], isCumSum = false): untyped =
  ## Computes the EDF of binned data
  var dataCdf = data
  if not isCumSum:
    seqmath.cumsum(dataCdf)
  let integral = dataCdf[^1]
  ## XXX: why min?
  let baseline = min(data) # 0.0 
  dataCdf.mapIt((it - baseline) / (integral - baseline))

proc sample*(cdf: seq[float], ys: seq[float]): float =
  let point = rand(1.0)
  let idx = cdf.lowerBound(point)
  if idx < cdf.len:
    result = ys[idx]
  else:
    result = Inf

proc expFn(x: float, λ: float): float =
  result = 1.0 / λ * exp(- x / λ)

proc sampleFrom*(fn: proc(x: float): float,
                low, high: float,
                num = 1000,
                samples = 1_000_000): seq[float] =
  ## Note: it may be useful to hand a closure with wrapped arguments!
  let xs = linspace(low, high, num)
  let ys = xs.mapIt( fn(it) ) 
  
  # now sample 100,000 points
  let cdf = ys.toEdf()
  result = toSeq(0 ..< samples).mapIt(sample(cdf, xs))

when isMainModule:
  ## Mini test: Compare with plot output from /tmp/test_sample.nim!
  let λ = 2.0
  
  let fnSample = (proc(x: float): float =
    result = expFn(x, λ)
  )
  
  let ySampled = sampleFrom(fnSample, 0.0, 3.0)
  let ySampled2 = sampleFrom(fnSample, 0.0, 10.0)
  proc toHisto(xs: seq[float]): DataFrame =
    const binSize = 0.1
    let binNum = ((xs.max - xs.min) / binSize).round.int      
    let (hist, bins) = histogram(xs, binNum)
    let maxH = hist.max
    result = toDf({"x" : bins[0 ..< ^2], "y" : hist.mapIt(it / maxH)})
  let dfC = bind_rows([("1", ySampled.toHisto()), ("2", ySampled2.toHisto())], "val")
  ggplot(dfC, aes("x", "y", fill = "val")) +
    #geom_histogram(bins = 100, density = true, alpha = 0.5, hdKind = hdOutline, fillColor = "red") +
    geom_histogram(bins = 100, alpha = 0.5, hdKind = hdOutline, stat = "identity", position = "identity") +
    ggsave("/t/test_sample_from.pdf")

import std / [random, sequtils, algorithm, strformat]
import seqmath, ggplotnim, unchained

import /tmp/sampling_helper

proc expFn(x: float, λ: float): float =
  result = 1.0 / λ * exp(- x / λ)

proc main(λ: float) =
  let σT = 640.0 # μm·√cm
  let fnSample = (proc(x: float): float =
    result = expFn(x, λ)
  )
  
  proc rmsTrans(x: float): float =
    let toDrift = (3.0 - x)
    result = sqrt(toDrift) * σT

  # sample from our exponential distribution describing absorption  
  let ySampled = sampleFrom(fnSample, 0.0, 3.0)
  # now compute the long and trans RMS for each
  let yRmsTrans = ySampled.mapIt(rmsTrans(it))
  ggplot(toDf(yRmsTrans), aes("yRmsTrans")) +
    geom_histogram(bins = 100, density = true, alpha = 0.5, hdKind = hdOutline, fillColor = "red") +
      ggsave(&"/t/sample_transverse_rms_{λ}_cm_absorption_length.pdf")
  
  #let sampleTransFn = (proc(x: float): float =
  #  result = gaus(x = x, mean = σT, 
when isMainModule:
  import cligen
  dispatch main

basti at voidRipper in /t λ ./simulate_rms_transverse_simple --λ 3.0 # below Ar absorption edge
basti at voidRipper in /t λ ./simulate_rms_transverse_simple --λ 2.2 # 5.9 keV
basti at voidRipper in /t λ ./simulate_rms_transverse_simple --λ 0.5 # 3.x keV above Ar absorption edge

cd /tmp/
mkdir RmsTransversePlots && cd RmsTransversePlots

plotData --h5file ~/CastData/data/CalibrationRuns2017_Reco.h5 \
 --runType rtCalibration \
 --chips 3 \
 --config ~/CastData/ExternCode/TimepixAnalysis/Plotting/karaPlot/config.toml \
 --ingrid \
 --cuts '("rmsTransverse", 0.1, 1.5)' \
 --cuts '("energyFromCharge", 2.5, 3.2)' \
 --applyAllCuts \
 --region crSilver

plotData --h5file ~/CastData/data/CalibrationRuns2017_Reco.h5 \
 --runType rtCalibration \
 --chips 3 \
 --config ~/CastData/ExternCode/TimepixAnalysis/Plotting/karaPlot/config.toml \
 --ingrid \
 --cuts '("rmsTransverse", 0.1, 1.5)' \
 --cuts '("energyFromCharge", 5.5, 6.5)' \
 --applyAllCuts \
 --region crSilver

cdl_spectrum_creation -i ~/CastData/data/CDL_2019/CDL_2019_Reco.h5 --dumpAccurate --hideNloptFit

import std / [random, sequtils, algorithm, strformat]
import seqmath, ggplotnim, unchained

import /tmp/sampling_helper

proc expFn(x: float, λ: float): float =
  result = 1.0 / λ * exp(- x / λ)

import random / mersenne
import alea / [core, rng, gauss]

proc main(E = 5.9, λ = 0.0) =
  ## Introduce sampling of a gaussian around σT with something like this
  ## which is ~150 = 1σ for a √3cm drift (seen in C-EPIC 0.6 kV CDL line
  ## rmsTransverse data)
  ## Note: another number we have for ΔσT is of course the simulation error
  ## on σT, but I suspect that's not a good idea (also it's large, but still
  ## much smaller than this).
  let ΔσT = 86.0 #  / 2.0

  ## XXX: Implement calculation of absorption length from `xrayAttenuation`
  # let dfAbs = 

  ## XXX: Implement extraction of diffusion values from data:
  let dfGas = readCsv("/home/basti/org/resources/ar_iso_97_7_2_3_septemboard_cast_different_temps.csv")

  let σT = 640.0 # μm/√cm
  let fnSample = (proc(x: float): float =
    result = expFn(x, λ)
  )


  var rnd = wrap(initMersenneTwister(1337))
  var gaus = gaussian(0.0, 1.0) # we will modify this gaussian for every draw!
  proc rmsTrans(x: float): float =
    let toDrift = (3.0 - x)
    # adjust the gaussian to Diffusion = σ_T · √(drift distance)
    # and width of Sigma = ΔσT · √(drift distance) (at 3 cm we want Δ of 150)
    gaus.mu = sqrt(toDrift) * σT
    gaus.sigma = ΔσT * sqrt(toDrift)
    #echo "DRAWING AROUND: ", gaus.mu, " WITH SIGMA: ", gaus.sigma
    result = rnd.sample(gaus)

  # sample from our exponential distribution describing absorption  
  let ySampled = sampleFrom(fnSample, 0.0, 3.0)
  # now compute the long and trans RMS for each
  let yRmsTrans = ySampled.mapIt(rmsTrans(it))
  let 
    GoldenMean = (sqrt(5.0) - 1.0) / 2.0  # Aesthetic ratio
    FigWidth = 1200.0                      # width in pixels
    FigHeight = FigWidth * GoldenMean     # height in pixels

  ggplot(toDf(yRmsTrans), aes("yRmsTrans")) +
    geom_histogram(bins = 100, density = true, alpha = 0.5, hdKind = hdOutline, fillColor = "red") +
    ggsave(&"/t/sample_gauss_transverse_rms_{λ}_cm_absorption_length.pdf", width = FigWidth, height = FigHeight)
  
when isMainModule:
  import cligen
  dispatch main

basti at voidRipper in /t λ ./simulate_rms_transverse_gauss --λ 3.0
basti at voidRipper in /t λ ./simulate_rms_transverse_gauss --λ 2.2
basti at voidRipper in /t λ ./simulate_rms_transverse_gauss --λ 2.0
basti at voidRipper in /t λ ./simulate_rms_transverse_gauss --λ 1.0
basti at voidRipper in /t λ ./simulate_rms_transverse_gauss --λ 0.5
basti at voidRipper in /t λ ./simulate_rms_transverse_gauss --λ 0.1

plotData --h5file ~/CastData/data/CDL_2019/CDL_2019_Reco.h5 \
 --runType rtCalibration \
 --chips 3 \
 --config ~/CastData/ExternCode/TimepixAnalysis/Plotting/karaPlot/config.toml \
 --ingrid \
 --cuts '("rmsTransverse", 0.1, 1.5)' \
 --applyAllCuts \
 --runs 342 --runs 343 \ 
 --region crSilver

plotData --h5file ~/CastData/data/CDL_2019/CDL_2019_Reco.h5 \
 --runType rtCalibration \
 --chips 3 \
 --config  ~/CastData/ExternCode/TimepixAnalysis/Plotting/karaPlot/config.toml \
 --ingrid \
 --cuts '("rmsTransverse", 0.1, 1.5)' \
 --applyAllCuts \
 --runs 328 --runs 329 --runs 351 \
 --region crSilver

plotData --h5file ~/CastData/data/CDL_2019/CDL_2019_Reco.h5 \
 --runType rtCalibration \
 --chips 3 \
 --config  ~/CastData/ExternCode/TimepixAnalysis/Plotting/karaPlot/config.toml \
 --ingrid \
 --cuts '("rmsTransverse", 0.1, 1.5)' \
 --applyAllCuts \
 --runs 328 --runs 329 --runs 351 \
 --region crSilver