Skip to content

analysis

get_reach(filename, threshold_meV=1.0, exposure_kg_yr=1.0, n_cut=3.0, model=None, time=0)

Computes the projected reach: the 95% C.L. constraint (3 events, no background) on \(ar{sigma}_n\) or \(ar{sigma}_e\) for a given model, in units of cm2 and normalizing to the appropriate reference cross section.

Parameters:

Name Type Description Default
filename str

The path to the HDF5 file containing the data.

 required
threshold_meV float

The threshold in meV. Defaults to 1.0.

 1.0
exposure_kg_yr float

The exposure in kg.yr. Defaults to 1.0.

 1.0
n_cut float

The number of events. Defaults to 3.0.

 3.0
model str

The model to use. Defaults to None (find it in file)

 None
time float | str

The time to use. Defaults to 0.

 0

Returns:

Type Description
array

np.array: the cross section.
Source code in darkmagic/analysis.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
def get_reach(
    filename: str,
    threshold_meV: float = 1.0,
    exposure_kg_yr: float = 1.0,
    n_cut: float = 3.0,
    model: str | None = None,
    time: float | str = 0,
) -> np.array:
    """
    Computes the projected reach: the 95% C.L. constraint (3 events, no background) on $\bar{sigma}_n$ or $\bar{sigma}_e$ for a given model, in units of cm2 and normalizing to the appropriate reference cross section.

    Args:
        filename (str): The path to the HDF5 file containing the data.
        threshold_meV (float, optional): The threshold in meV. Defaults to 1.0.
        exposure_kg_yr (float, optional): The exposure in kg.yr. Defaults to 1.0.
        n_cut (float, optional): The number of events. Defaults to 3.0.
        model (str, optional): The model to use. Defaults to None (find it in file)
        time (float | str, optional): The time to use. Defaults to 0.

    Returns:
        np.array: the cross section.
    """

    with h5py.File(filename, "r") as f:
        energy_bin_width = f["numerics"]["energy_bin_width"][...]
        threshold = f["particle_physics"]["threshold"][...]
        m_chi = np.array(f["particle_physics"]["dm_properties"]["mass_list"])
        raw_binned_rate = 1e-100 + np.array(  # 1e-100 to avoid division by zero
            [
                np.array(f["data"]["diff_rate"][str(time)][str(m)])
                for m in range(len(m_chi))
            ]
        )
        # TODO: add check to see if model is there

    # Vanilla PhonoDark, have a threshold of 1 meV by default
    # We need to account for that in the binning
    bin_num_cut = int(threshold_meV * 1e-3 / energy_bin_width) - int(
        threshold / energy_bin_width
    )

    cs_constraint = n_cut / np.sum(raw_binned_rate[:, bin_num_cut:], axis=1)
    cs_constraint /= (
        exposure_kg_yr
        * const.kg_yr
        * const.cm2
        * BUILT_IN_MODELS[model].ref_cross_sect(m_chi)
    )

    return cs_constraint