Add method to calculate the exposure (#75)

src/LegendDataManagement.jl

@@ -50,6 +50,7 @@ include("evt_functions.jl")
 include("lprops.jl")
 include("data_io.jl")
 include("active_volume.jl")
+include("exposure.jl")
 include("utils/utils.jl")
 
 end # module

src/exposure.jl

@@ -0,0 +1,112 @@
+# This file is a part of LegendDataManagement.jl, licensed under the MIT License (MIT).
+
+"""
+    get_exposure(data::LegendData, det::DetectorIdLike, period::DataPeriodLike, run::DataRunLike; kwargs...)
+    get_exposure(data::LegendData, det::DetectorIdLike, period::DataPeriodLike; kwargs...)
+    get_exposure(data::LegendData, det::DetectorIdLike, partition::DataPartitionLike; kwargs...)
+
+Calculates the exposure of a detector in a given run/period/partition.
+
+# Arguments
+- `data`: `LegendData` object with information on detector geometries and `runinfo` / `partitioninfo`
+- `det` Detector for which the exposure is calculated
+- `period`: `DataPeriod` for which the exposure is calculated
+- `run`: `DataRun` for which the exposure is calculated
+- `partition`: `DataPartition` for which the exposure is calculated
+
+# Keyword Arguments
+- `is_analysis_run`: If set to `true`, only the `runs` flagged as `is_analysis_phy_tun == true` are considered. Default is `true`.
+- `cat` `DataCategory` for which the exposure is calculated. Default is `:phy`.`
+
+# Returns
+- `exposure`: the exposure of the detector `det` for the time given.
+
+
+# Example
+```julia
+l200 = LegendData(:l200)
+get_exposure(l200, :V00050A, DataPeriod(3), DataRun(0))
+get_exposure(l200, :V00050A, DataPeriod(3))
+get_exposure(l200, :V00050A, DataPartition(1))
+````
+
+"""
+function get_exposure(data::LegendData, det::DetectorIdLike, period::DataPeriodLike, run::DataRunLike; is_analysis_run::Bool=true, cat::DataCategoryLike=:phy)
+    rinfo = runinfo(data, period, run)
+    _get_exposure(data, det, Table([rinfo]), is_analysis_run, cat)
+end
+
+function get_exposure(data::LegendData, det::DetectorIdLike, period::DataPeriod; is_analysis_run::Bool=true, cat::DataCategoryLike=:phy)
+    rinfo = runinfo(data, period)
+    _get_exposure(data, det, rinfo, is_analysis_run, cat)
+end
+
+function get_exposure(data::LegendData, det::DetectorIdLike, part::DataPartition; is_analysis_run::Bool=true, cat::DataCategoryLike=:phy)
+    part_dict = partitioninfo(data, det)
+    if haskey(part_dict, part)
+        rinfo = partitioninfo(data, det, part)
+        return _get_exposure(data, det, rinfo, is_analysis_run, cat)
+    end
+
+    #default if partition does not exist
+    return 0.0u"kg*yr"
+end
+
+function get_exposure(data::LegendData, det::DetectorIdLike, sel::Union{AbstractString, Symbol}; kwargs...)
+    selectors = (DataPartition, DataPeriod)
+    for SEL in selectors
+        if _can_convert_to(SEL, sel)
+             return _get_exposure(data, det, SEL(sel); kwargs...)
+        end
+    end
+    throw(ArgumentError("The selector $(sel) cannot be converted to type: $(selectors)"))
+end
+
+
+### TODO: determine livetimes from data files instead of metadata
+function _get_exposure(data::LegendData, det::DetectorIdLike, rinfo::Table, is_analysis_run::Bool=true, cat::DataCategoryLike=:phy)
+
+    # check that the DataCategory is valid
+    if !(_can_convert_to(DataCategory, cat) && hasproperty(rinfo, DataCategory(cat).label))
+        throw(ArgumentError("Data category `$(cat)`` is invalid"))
+    end
+    cat_label::Symbol = DataCategory(cat).label
+
+    # determine livetime
+    rinfo_cat = getproperty(rinfo, cat_label)
+    livetimes = getproperty.(rinfo_cat, :livetime)
+
+    # if is_analysis_run == true: 
+    # check that the analysis flag is valid and apply it
+    analysis_flag = Symbol("is_analysis_$(cat_label)_run")
+    if is_analysis_run 
+        if !hasproperty(rinfo, analysis_flag)
+            throw(ArgumentError("No column `$(analysis_flag)` found. Please set `is_analysis_run = false` in `get_exposure`"))
+        end
+        livetimes = livetimes .* getproperty(rinfo, analysis_flag)
+    end
+    # sum up all livetimes (excluding NaN values)
+    livetime = !isempty(livetimes) ? sum((livetimes .* .!isnan.(livetimes))) : 0.0u"s"
+
+    # determine the mass of 76Ge
+    filekeys = getproperty.(rinfo_cat, :startkey)
+    mass = if !iszero(livetime) && !isempty(filekeys)
+        # read in the channelinfo
+        filekey = first(filekeys)
+        _chinfo = channelinfo(data, filekey, det, extended = true, verbose = false)
+        chinfo = if !all(x -> hasproperty(_chinfo, x), (:enrichment, :mass, :active_volume, :total_volume))
+            empty!(_cached_channelinfo)
+            channelinfo(data, filekey, det, extended = true, verbose = false)
+        else
+            _chinfo
+        end
+        # chinfo.active_volume == chinfo.total_volume && @warn "No FCCD value given for detector $(det)"
+        chinfo.mass * chinfo.enrichment * chinfo.active_volume / chinfo.total_volume
+    else
+        0.0u"kg"
+    end
+
+    return uconvert(u"kg*yr", livetime * mass)
+end
+
+export get_exposure

src/legend_data.jl

@@ -245,7 +245,7 @@ const _cached_channelinfo = LRU{Tuple{UInt, AnyValiditySelection}, StructVector}
 Get all channel information for the given [`LegendData`](@ref) and
 [`ValiditySelection`](@ref).
 """
-function channelinfo(data::LegendData, sel::AnyValiditySelection; system::Symbol = :all, only_processable::Bool = false, extended::Bool = false)
+function channelinfo(data::LegendData, sel::AnyValiditySelection; system::Symbol = :all, only_processable::Bool = false, extended::Bool = false, verbose::Bool = true)
     key = (objectid(data), sel)
     chinfo = get!(_cached_channelinfo, key) do
         chmap = data.metadata(sel).hardware.configuration.channelmaps
@@ -320,7 +320,7 @@ function channelinfo(data::LegendData, sel::AnyValiditySelection; system::Symbol
                                    isa(fccds, PropDicts.MissingProperty) || 
                                    isa(fccds[first(keys(fccds))].value, PropDicts.MissingProperty)
 
-                    haskey(diodmap, k) && @warn "No FCCD value given for detector $(detector)"
+                    verbose && haskey(diodmap, k) && @warn "No FCCD value given for detector $(detector)"
                     0.0
                 else 
                     fccds[first(keys(fccds))].value

test/Project.toml

@@ -17,3 +17,4 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 Documenter = "1"
+LegendTestData = "0.2.9"

test/runtests.jl

@@ -15,6 +15,7 @@ Test.@testset "Package LegendDataManagement" begin
     include("test_lpy_expressions.jl")
     include("test_dataprod_config.jl")
     include("test_lprops.jl")
+    include("test_exposure.jl")
     include("test_ext_ssd.jl")
     include("test_ext_plots.jl")
     include("test_ext_legendhdf5io.jl")

test/test_exposure.jl

@@ -0,0 +1,31 @@
+# This file is a part of LegendDataManagement.jl, licensed under the MIT License (MIT).
+
+using Test
+using LegendDataManagement
+using Unitful
+
+include("testing_utils.jl")
+
+l200 = LegendData(:l200)
+@testset "Exposure" begin
+    for det in (:V99000A, :B99000A)
+        @testset "$(det)" begin
+            @testset "Period exposure" begin
+                period = DataPeriod(2)
+                rinfo = runinfo(l200, period)
+                period_exposure = get_exposure(l200, det, period)
+                @test period_exposure isa Quantity
+                @test dimension(period_exposure) == dimension(u"kg*yr")
+                @test period_exposure ≈ sum(map(r -> get_exposure(l200, det, period, r), rinfo.run))
+            end
+            @testset "Partition exposure" begin
+                part = DataPartition(1)
+                part_exposure = get_exposure(l200, det, part)
+                partinfo = partitioninfo(l200, det, part)
+                @test part_exposure isa Quantity
+                @test dimension(part_exposure) == dimension(u"kg*yr")
+                @test part_exposure ≈ sum(map(p -> get_exposure(l200, det, p.period, p.run), partinfo))  
+            end
+        end
+    end
+end