import shell, strutils

proc main(path: string, run: int) =
  # parse data
  let outfile = "/t/xray_finger_$#.h5" % $run
  let recoOut = "/t/reco_xray_finger_$#.h5" % $run
  
  shell:
    raw_data_manipulation -p ($path) "--runType xray --out " ($outfile)
  shell:
    reconstruction -i ($outfile) "--out " ($recoOut)
  
when isMainModule:
  import cligen
  dispatch main

import nimhdf5, ggplotnim, options
import ingrid / tos_helpers
import std / [strutils, tables]

proc main(run: int, switchAxes: bool = false, useTeX = false) =
  let file = "/t/reco_xray_finger_$#.h5" % $run
  
  #proc readClusters(h5f: H5File): (seq[float], seq[float]) =
  var h5f = H5open(file, "r")
  
  # compute counts based on number of each pixel hit
  proc toIdx(x: float): int = (x / 14.0 * 256.0).round.int.clamp(0, 255)
  var ctab = initCountTable[(int, int)]() 
  
  var df = readRunDsets(h5f, run = run,
                        chipDsets = some((
                        chip: 3, dsets: @["centerX", "centerY"])))
    .mutate(f{"xidx" ~ toIdx(idx("centerX"))},
            f{"yidx" ~ toIdx(idx("centerY"))})
  let xidx = df["xidx", int]
  let yidx = df["yidx", int]
  forEach x in xidx, y in yidx:
    inc cTab, (x, y)
  df = df.mutate(f{int: "count" ~ cTab[(`xidx`, `yidx`)]})
  let centerX = df["centerX", float].mean
  let centerY = df["centerY", float].mean
  discard h5f.close()
 
  echo "Center position of the cluster is at: (x, y) = (", centerX, ", ", centerY, ")"
  ## NOTE: Exchanging the axes for X and Y is equivalent to a 90° clockwise rotation for our data
  ## because the centerX values are inverted `(256 - x), applyPitchConversion`. 
  ## The real rotation of the Septemboard detector at CAST seen from the telescope onto the
  ## detector is precisely 90° clockwise. 
  let x = if switchAxes: "centerY" else: "centerX"
  let y = if switchAxes: "centerX" else: "centerY"
  let cX = if switchAxes: centerY else: centerX
  let cY = if switchAxes: centerX else: centerY
  ggplot(df, aes(x, y, color = "count")) +
    geom_point(size = 0.75) +
    geom_point(data = newDataFrame(), aes = aes(x = cX, y = cY),
               color = "red", marker = mkRotCross) + 
    scale_color_continuous() +
    ggtitle("X-ray finger clusters of run $#" % $run) +
    xlab(r"x [mm]") + ylab(r"y [mm]") + 
    xlim(0.0, 14.0) + ylim(0.0, 14.0) +
    margin(right = 3.5) + 
    #theme_scale(1.0, family = "serif") +
    coord_fixed(1.0) + 
    themeLatex(fWidth = 0.5, width = 600, baseTheme = sideBySide) +
    legendPosition(0.83, 0.0) + 
    ggsave("/home/basti/phd/Figs/CAST_Alignment/xray_finger_centers_run_$#.pdf" % $run,
           useTeX = useTeX, standalone = useTeX, dataAsBitmap = true)
           #useTeX = true, standalone = true)

when isMainModule:
  import cligen
  dispatch main

./code/xray_finger_data_parsing -p ~/CastData/data/XrayFingerRuns/Run_21_170713-11-03 --run 21
./code/xray_finger_data_parsing -p ~/CastData/data/XrayFingerRuns/Run_189_180420-09-53 --run 189

./code/xray_finger_center_plot -r 21 --switchAxes 
./code/xray_finger_center_plot -r 189 --switchAxes

plotBackgroundRate \
    \t\reco_xray_finger_189.h5 \
    --names "X-ray finger" \
    --title "X-ray finger run 189 spectrum" \
    --centerChip 3 \
    --region crGold \
    --energyDset energyFromCharge \
    --outfile xray_finger_spectrum_189.pdf \
    --outpath ~/phd/Figs/XrayFinger/ \
    --useTeX \
    --quiet

import tables

type
  # defines the TubesMap datatype, which is a combined object to
  # store the different parts of the tubing each sequences of tuples
  TubesMap* = object
    static_tubes* : seq[tuple[diameter: float, length: float]]
    flexible_tubes* : seq[tuple[diameter: float, length: float]]
    t_pieces* : seq[tuple[diameter: float, length_long: float, length_short: float]]
    crosses* : seq[tuple[diameter: float, length: float]]

proc getVacuumTubing*(): TubesMap =
  # this function returns the data (originally written in calc_vacuum_volume.org
  # as a set of hash maps as a "TubesMap" datatype
  let st_tubing = @[(63.0, 10.0),
                    (63.0, 51.0),
                    (63.0, 21.5),
                    (25.0, 33.7),
                    (63.0, 20.0),
                    (63.0, 50.0),
                    (40.0, 15.5),
                    (16.0, 13.0),
                    (40.0, 10.0)]

  let fl_tubing = @[(16.0,  25.0),
                    (16.0,  25.0 ),
                    (16.0,  25.0 ),
                    (16.0,  25.0 ),
                    (16.0,  40.0 ),
                    (25.0,  90.0 ),
                    (25.0,  80.0 ),
                    (40.0,  50.0 ),
                    (16.0, 150.0 ),
                    (40.0,  80.0 ),
                    (40.0,  80.0)]

  let t_pieces = @[(40.0, 18.0, 21.0),
                   (16.0, 7.0, 4.5),
                   (40.0, 10.0, 10.0)]

  let crosses = @[(16.0, 10.0),
                  (40.0, 14.0),
                  (40.0, 14.0),
                  (40.0, 14.0)]
                    
  let t = TubesMap(static_tubes: st_tubing, flexible_tubes: fl_tubing, t_pieces: t_pieces, crosses: crosses)
  echo "Vacuum tubing is as follows:"
  echo t
  return t

import math
import tubing
import sequtils, future
import typeinfo

# This script contains a calculation for the total volume of the
# currently in use vacuum system at CAST (behind and including LLNL
# telescope)

proc cylinder_volume(diameter, length: float): float =
  # this proc calculates the volume of a cylinder, given a
  # diameter and a length both in cm
  result = PI * pow(diameter / 2.0, 2) * length
  
proc t_piece_volume(diameter, length_long, length_short: float): float =
  # this proc calculates the volume of a T shaped vacuum piece, using
  # the cylinder volume proc
  # inputs:
  # diameter: diameter of the tubing in cm
  # length_long: length of the long axis of the tubing
  # length_short: length of the short axis of the tubing
  result = cylinder_volume(diameter, length_long) + cylinder_volume(diameter, length_short - diameter)

proc cross_piece_volume(diameter, length: float): float =
  # this proc calculates the volume of a cross shaped vacuum piece, using
  # the cylinder volume proc
  # inputs:
  # diameter: diameter of the tubing in cm
  # length: length of one axis of the tubing
  result = 2 * cylinder_volume(diameter, length) - pow(diameter, 3)

proc calcTotalVacuumVolume(t: TubesMap): float =
  # function which calculates the total vacuum volume, using
  # the rough measurements of the length and diameters of all the
  # piping
  # the TubesMap consists of:
  # static_tubes : seq[tuple[diameter: float, length: float]]
  # flexible_tubes : seq[tuple[diameter: float, length: float]]
  # t_pieces : seq[tuple[diameter: float, length_long: float, length_short: float]]
  # crosses : seq[tuple[diameter: float, length: float]]
  # define variables to store static volume etc

  # calc volume of static tubing
  let static_vol = sum(map(
    t.static_tubes, (b: tuple[diameter, length: float]) ->
    float =>
    cylinder_volume(b.diameter / 10, b.length)))
  let flexible_vol = sum(map(
    t.flexible_tubes, (b: tuple[diameter, length: float]) -> 
    float =>
    cylinder_volume(b.diameter / 10, b.length)))
  let t_vol = sum(map(
    t.t_pieces, (b: tuple[diameter, length_long, length_short: float]) ->
    float =>
    t_piece_volume(b.diameter / 10, b.length_long, b.length_short)))
  let crosses_vol = sum(map(
    t.crosses, (b: tuple[diameter, length: float]) ->
    float =>
    cross_piece_volume(b.diameter / 10, b.length)))

  result = static_vol + flexible_vol + t_vol + crosses_vol

proc calcFlowRate(d, p, mu, x: float): float =
  # this function calculates the flow rate following the Poiseuille Equation
  # for a non-ideal gas under laminar flow.
  # inputs:
  # d: diameter of the tube in m
  # p: pressure difference between both ends of the tube in Pa
  # mu: dynamic viscosity of the medium
  # x: length of the tube
  # note: get viscosity e.g. from https://www.lmnoeng.com/Flow/GasViscosity.php
  # returns the flow rate in m^3 / s
  result = PI * pow(d, 4) * p / (128 * mu * x)

proc calcGasAmount(p, V, T: float): float =
  # this function calculates the amount of gas in moles follinwg
  # the ideal gas equation p V = n R T for a given pressure, volume
  # and temperature
  let R = 8.31446
  result = p * V / (R * T)

proc calcVolumeFromMol(p, n, T: float): float =
  # this function calculates the volume in m^3 follinwg
  # the ideal gas equation p V = n R T for a given pressure, amount in mol
  # and temperature
  let R = 8.31446
  result = n * R * T / p
    
proc main() =

  # TODO: checke whether diameter of 63mm for telescope is a reasonable
  # number!
  let t = getVacuumTubing()
  # first of all we need to calculate the total volume of the vacuum
  let volume = calcTotalVacuumVolume(t)
  echo volume

  # now calcualte flow rate through pipe
  let
    # 3 mm diameter
    d = 3e-3
    # 6 bar pressure diff
    p = 6.0e5
    # viscosity of air
    mu = 1.8369247e-4
    # ~2m of tubing
    x = 2.0
    flow = calcFlowRate(d, p, mu, x)

  echo(flow * 1e3, " l / s")

  # given the flow in liter, calc total gas inserted into the system
  let flow_l = flow * 1e3

  # detector volume in m^3
  let det_vol = cylinder_volume(12.0, 3.0) * 1e-6
  echo("Detector volume is : ", det_vol)
  # initial gas volume inside detector (1 bar is argon!), thus
  # only .5 bar
  let n_initial = calcGasAmount(0.5e5, det_vol, 293.15)
  # gas which came in after window ruptured
  let valve_open = 5.0
  # total volume in m^3
  let flow_vol = flow_l * 1e-3 * valve_open
  
  # since the flown volume is given for normal pressure and temp, calc
  # amount of gas
  let n_flow = calcGasAmount(1.0e5, flow_vol, 293.15)
  echo("Initial gas is : ", n_initial, " mol")
  echo("Gas from flow is : ", n_flow, " mol")
  let n_total = n_initial + n_flow
  echo("Total compressed air, which entered system : ", n_total)

  # calc volume corresponding to normal pressure
  let tot_vol_atm = calcVolumeFromMol(1e5, n_total, 293.15)
  echo("Total volume of air at normal pressure : ", tot_vol_atm * 1e3, " l")
  
when isMainModule:
  main()

# 1. open file
# 2. get file info
# 3. for each run get extended run info
# ?
import std / [times, strformat, strutils]
import nimhdf5, unchained
import ingrid / tos_helpers

type
  CastInformation = object
    totalDuration: Second
    activeDuration: Second
    activeFraction: float # ratio of active / total    
    numTrackings: int
    nonTrackingDuration: Second
    activeNonTrackingTime: Second
    trackingTime: Second
    activeTrackingTime: Second
    totalEvents: int # total number of recorded events
    onlyCenter: int  # events with activity only on center chip (> 3 hits)
    onlyOuter: int   # events with activity only on outer, but not center chip
    centerAndOuter: int # events with activity on center & any outer chip
    center: int      # events with activity on center (irrespective any other)
    anyActive: int # events with any active chip
    fractionWithCenter: float # fraction of events that have center chip activity
    fractionWithAny: float # fraction of events that have any activity
    # ... add mean of event durations?
    fadcReadouts: int
    fractionFadc: float # fraction of events having FADC readout
    scinti1NonTrivial: int # number of non trivial scinti triggers 0 < x < 4095
    scinti2NonTrivial: int # number of non trivial scinti triggers 0 < x < 4095
    scinti1Triggers: int # number of any scinti triggers != 0
    scinti2Triggers: int # number of any scinti triggers != 0
    fractionScinti1: float # fraction of events with any scinti 1 activity    
    fractionScinti2: float # fraction of events with any scinti 2 activity

proc fieldToStr(s: string): string =
  case s
  of "totalDuration":         result = "total duration"
  of "activeDuration":        result = "active duration"
  of "activeFraction":        result = "active fraction"
  of "numTrackings":          result = "# trackings"
  of "nonTrackingDuration":   result = "non tracking time"
  of "activeNonTrackingTime": result = "active non tracking time"
  of "trackingTime":          result = "tracking time"
  of "activeTrackingTime":    result = "active tracking time"
  of "totalEvents":           result = "total # events"
  of "center":                result = "center chip"
  of "onlyCenter":            result = "only center chip"
  of "onlyOuter":             result = "only any outer chip"
  of "centerAndOuter":        result = "center + outer"
  of "anyActive":             result = "any chip"
  of "fractionWithCenter":    result = "fraction with center"
  of "fractionWithAny":       result = "fraction with any"
  of "fadcReadouts":          result = "with fadc readouts"
  of "fractionFadc":          result = "fraction with FADC"
  of "scinti1NonTrivial":     result = "with SiPM trigger <4095"
  of "scinti2NonTrivial":     result = "with veto scinti trigger <4095"
  of "scinti1Triggers":       result = "with any SiPM trigger"
  of "scinti2Triggers":       result = "with any veto scinti trigger"
  of "fractionScinti1":       result = "fraction with any SiPM"
  of "fractionScinti2":       result = "fraction with any veto scinti"  

proc `$`(castInfo: CastInformation): string =
  result.add &"Total duration: {pretty(castInfo.totalDuration.to(Hour), 4, true)}\n"
  result.add &"Active duration: {pretty(castInfo.activeDuration.to(Hour), 4, true)}\n"
  result.add &"Active fraction: {castInfo.activeFraction}\n"
  result.add &"Number of trackings: {castInfo.numTrackings}\n"
  result.add &"Non-tracking time: {pretty(castInfo.nonTrackingDuration.to(Hour), 4, true)}\n"
  result.add &"Active non-tracking time: {pretty(castInfo.activeNonTrackingTime.to(Hour), 4, true)}\n"
  result.add &"Tracking time: {pretty(castInfo.trackingTime.to(Hour), 4, true)}\n"
  result.add &"Active tracking time: {pretty(castInfo.activeTrackingTime.to(Hour), 4, true)}\n"
  result.add &"Number of total events:   {castInfo.totalEvents}\n"
  result.add &"Number of events without center: {castInfo.onlyOuter}\n"
  result.add &"\t| {(castInfo.onlyOuter.float / castInfo.totalEvents.float) * 100.0} %\n"
  result.add &"Number of events only center: {castInfo.onlyCenter}\n"
  result.add &"\t| {(castInfo.onlyCenter.float / castInfo.totalEvents.float) * 100.0} %\n"
  result.add &"Number of events with center activity and outer: {castInfo.centerAndOuter}\n"
  result.add &"\t| {(castInfo.centerAndOuter.float / castInfo.totalEvents.float) * 100.0} %\n"
  result.add &"Number of events any hit events: {castInfo.anyActive}\n"
  result.add &"\t| {(castInfo.anyActive.float / castInfo.totalEvents.float) * 100.0} %\n"

proc countEvents(df: DataFrame): int =
  for (tup, subdf) in groups(df.group_by("runNumber")):
    inc result, subDf["eventNumber", int].max

proc contains[T](t: Tensor[T], x: T): bool =
  for i in 0 ..< t.size:
    if x == t[i]:
      return true

proc countChipActivity(castInfo: var CastInformation, df: DataFrame) =
  for (tup, subDf) in groups(df.group_by(["eventNumber", "runNumber"])):
    let chips = subDf["chip"].unique.toTensor(int)
    if 3 in chips:
      inc castInfo.center
    # start new if 
    if 3 notin chips: 
      inc castInfo.onlyOuter
    elif [3].toTensor == chips:
      inc castInfo.onlyCenter
    elif 3 in chips and chips.len > 1:
      inc castInfo.centerAndOuter
    inc castInfo.anyActive

proc processFile(fname: string): CastInformation = # extend to both calib & both background
  let h5f = H5open(fname, "r")
  let fileInfo = getFileInfo(h5f)
  var castInfo: CastInformation
  for run in fileInfo.runs:
    let runInfo = getExtendedRunInfo(h5f, run, fileInfo.runType)
    castInfo.totalDuration  += runInfo.timeInfo.t_length.inSeconds().Second
    castInfo.activeDuration += runInfo.activeTime.inSeconds.Second
    castInfo.nonTrackingDuration += runInfo.nonTrackingDuration.inSeconds.Second
    castInfo.activeNonTrackingTime += runInfo.activeNonTrackingTime.inSeconds.Second    
    castInfo.numTrackings   += runInfo.trackings.len
    castInfo.trackingTime   += runInfo.trackingDuration.inSeconds.Second
    castInfo.activeTrackingTime += runInfo.activeTrackingTime.inSeconds.Second
    # read the data of all chips & FADC
    const names = ["eventNumber", "fadcReadout", "szint1ClockInt", "szint2ClockInt"]
    let dfNoChips = h5f.readRunDsets(run, commonDsets = names)
    let dfChips = h5f.readRunDsetsAllChips(run, fileInfo.chips,
                                           dsets = @[]) # don't need additional dsets
    castInfo.totalEvents       += dfNoChips.countEvents()
    castInfo.countChipActivity(dfChips)
    castInfo.fadcReadouts      += dfNoChips.filter(f{`fadcReadout` == 1}).len
    castInfo.scinti1Triggers   += dfNoChips.filter(f{`szint1ClockInt` != 0}).len
    castInfo.scinti2Triggers   += dfNoChips.filter(f{`szint2ClockInt` != 0}).len
    castInfo.scinti1NonTrivial += dfNoChips.filter(f{`szint1ClockInt` != 0 and `szint1ClockInt` < 4095}).len
    castInfo.scinti2NonTrivial += dfNoChips.filter(f{`szint2ClockInt` != 0 and `szint2ClockInt` < 4095}).len
      

# compute at the end as we need total information about fraction of total / active
  template fraction(arg, by: untyped): untyped = (castInfo.arg / castInfo.by) * 100.0
  castInfo.activeFraction = fraction(activeDuration, totalDuration)
  # fractions
  castInfo.fractionWithCenter = fraction(center         , totalEvents)
  castInfo.fractionWithAny    = fraction(anyActive      , totalEvents)
  castInfo.fractionFadc       = fraction(fadcReadouts   , totalEvents)
  castInfo.fractionScinti1    = fraction(scinti1Triggers, totalEvents)
  castInfo.fractionScinti2    = fraction(scinti2Triggers, totalEvents)
  echo castInfo
  result = castInfo

proc toTable(castInfos: Table[(string,string), CastInformation]): string =
  ## Turns the input into an Org table
  # | Field | Back Run-2 | Back Run-3 | Calib Run-2 | Calib Run-3 |
  # |-
  # ...
  # turn the input into a DF, then `toOrgTable` it
  proc toColName(tup: (string, string)): string =
    result = tup[1] & " "
    if "2017" in tup[0]:
      result.add "Run-2"
    else:
      result.add "Run-3"
      
  var df = newDataFrame()
  for k, v in pairs(castInfos):
    var fields = newSeq[string]()
    var vals = newSeq[string]()
    for field, val in fieldPairs(v):
      fields.add field.fieldToStr()
      when typeof(val) is Second:
        vals.add pretty(val.to(Hour), precision = 2, short = true,
                        format = ffDecimal)
      elif typeof(val) is float:
        vals.add $(val.formatFloat(precision = 4)) & " %"
      else:
        vals.add "\\num{" & $val & "}"
    let colName = k.toColName()
    let dfLoc = toDf({"Field" : fields, colName : vals})
    if df.len == 0:
      df = dfLoc
    else:
      df[colName] = dfLoc[colName]

  df = df.select(["Field", "calib Run-2", "calib Run-3", "back Run-2", "back Run-3"])
  echo df.toOrgTable(emphStrNumber = false)

proc main(background: seq[string], calibration: seq[string]) =

  var tab = initTable[(string, string), CastInformation]()
  for b in background:
    echo "--------------- Processing: ", b, " ---------------"
    tab[(b, "back")] = processFile(b)
    
  for c in calibration:
    echo "--------------- Processing: ", c, " ---------------"    
    tab[(c, "calib")] = processFile(c)
  echo tab
  echo tab.toTable()

when isMainModule:
  import cligen
  dispatch main

code/cast_run_information \
    ~/CastData/data/DataRuns2017_Reco.h5 \
    ~/CastData/data/DataRuns2018_Reco.h5 \
    -c ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    -c ~/CastData/data/CalibrationRuns2018_Reco.h5