Skip to content

Commit

Permalink
Merge pull request #8 from vkucera/cross-section
Browse files Browse the repository at this point in the history 
Fix the script utilities
  • Loading branch information
@DelloStritto
DelloStritto authored Apr 10, 2024
2 parents 3a3aae8 + e5f95ae commit e83093f
Show file tree
Hide file tree
Showing 2 changed files with 49 additions and 177 deletions.
154 changes: 31 additions & 123 deletions PWGHF/D2H/Macros/hf_analysis_utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,14 @@
author: Fabrizio Grosa <[email protected]>, CERN
"""

import numpy as np
import numpy as np # pylint: disable=import-error


def make_list(object) -> list:
"""
Returns the object as a list if it is not a list already.
"""
return object if isinstance(object, list) else list(object)


# pylint: disable=too-many-arguments
Expand Down Expand Up @@ -46,12 +53,7 @@ def compute_crosssection(
- crosssec_unc: cross-section statistical uncertainty
"""

crosssection = (
rawy
* frac
* sigma_mb
/ (2 * delta_pt * delta_y * eff_times_acc * n_events * br)
)
crosssection = rawy * frac * sigma_mb / (2 * delta_pt * delta_y * eff_times_acc * n_events * br)
if method_frac == "Nb":
crosssec_unc = rawy_unc / (rawy * frac) * crosssection
else:
Expand All @@ -68,7 +70,7 @@ def compute_fraction_fc(
cross_sec_fd,
raa_prompt=1.0,
raa_fd=1.0,
):
) -> "tuple[list[float], list[float]]":
"""
Method to get fraction of prompt / FD fraction with fc method
Expand All @@ -89,16 +91,13 @@ def compute_fraction_fc(
- frac_fd: list of fraction of non-prompt D (central, min, max)
"""

if not isinstance(cross_sec_prompt, list) and isinstance(cross_sec_prompt, float):
cross_sec_prompt = [cross_sec_prompt]
if not isinstance(cross_sec_fd, list) and isinstance(cross_sec_fd, float):
cross_sec_fd = [cross_sec_fd]
if not isinstance(raa_prompt, list) and isinstance(raa_prompt, float):
raa_prompt = [raa_prompt]
if not isinstance(raa_fd, list) and isinstance(raa_fd, float):
raa_fd = [raa_fd]
cross_sec_prompt = make_list(cross_sec_prompt)
cross_sec_fd = make_list(cross_sec_fd)
raa_prompt = make_list(raa_prompt)
raa_fd = make_list(raa_fd)

frac_prompt, frac_fd = [], []
frac_prompt: list[float] = []
frac_fd: list[float] = []
if acc_eff_prompt == 0:
frac_fd_cent = 1.0
frac_prompt_cent = 0.0
Expand All @@ -115,37 +114,11 @@ def compute_fraction_fc(
for i_sigma, (sigma_p, sigma_f) in enumerate(zip(cross_sec_prompt, cross_sec_fd)):
for i_raa, (raa_p, raa_f) in enumerate(zip(raa_prompt, raa_fd)):
if i_sigma == 0 and i_raa == 0:
frac_prompt_cent = 1.0 / (
1 + acc_eff_fd / acc_eff_prompt * sigma_f / sigma_p * raa_f / raa_p
)
frac_fd_cent = 1.0 / (
1 + acc_eff_prompt / acc_eff_fd * sigma_p / sigma_f * raa_p / raa_f
)
frac_prompt_cent = 1.0 / (1 + acc_eff_fd / acc_eff_prompt * sigma_f / sigma_p * raa_f / raa_p)
frac_fd_cent = 1.0 / (1 + acc_eff_prompt / acc_eff_fd * sigma_p / sigma_f * raa_p / raa_f)
else:
frac_prompt.append(
1.0
/ (
1
+ acc_eff_fd
/ acc_eff_prompt
* sigma_f
/ sigma_p
* raa_f
/ raa_p
)
)
frac_fd.append(
1.0
/ (
1
+ acc_eff_prompt
/ acc_eff_fd
* sigma_p
/ sigma_f
* raa_p
/ raa_f
)
)
frac_prompt.append(1.0 / (1 + acc_eff_fd / acc_eff_prompt * sigma_f / sigma_p * raa_f / raa_p))
frac_fd.append(1.0 / (1 + acc_eff_prompt / acc_eff_fd * sigma_p / sigma_f * raa_p / raa_f))

if frac_prompt and frac_fd:
frac_prompt.sort()
Expand All @@ -172,7 +145,7 @@ def compute_fraction_nb(
sigma_mb,
raa_ratio=1.0,
taa=1.0,
):
) -> "list[float]":
"""
Method to get fraction of prompt / FD fraction with Nb method
Expand All @@ -196,102 +169,39 @@ def compute_fraction_nb(
- frac: list of fraction of prompt (non-prompt) D (central, min, max)
"""

if not isinstance(crosssection, list) and isinstance(crosssection, float):
crosssection = [crosssection]

if not isinstance(raa_ratio, list) and isinstance(raa_ratio, float):
raa_ratio = [raa_ratio]
crosssection = make_list(crosssection)
raa_ratio = make_list(raa_ratio)

frac = []
frac: list[float] = []
for i_sigma, sigma in enumerate(crosssection):
for i_raa_ratio, raa_rat in enumerate(raa_ratio):
raa_other = 1.0
if i_sigma == 0 and i_raa_ratio == 0:
if raa_rat == 1.0 and taa == 1.0: # pp
frac_cent = (
1
- sigma
* delta_pt
* delta_y
* acc_eff_other
* br
* n_events
* 2
/ rawy
/ sigma_mb
)
frac_cent = 1 - sigma * delta_pt * delta_y * acc_eff_other * br * n_events * 2 / rawy / sigma_mb
else: # p-Pb or Pb-Pb: iterative evaluation of Raa needed
delta_raa = 1.0
while delta_raa > 1.0e-3:
raw_fd = (
taa
* raa_rat
* raa_other
* sigma
* delta_pt
* delta_y
* acc_eff_other
* br
* n_events
* 2
taa * raa_rat * raa_other * sigma * delta_pt * delta_y * acc_eff_other * br * n_events * 2
)
frac_cent = 1 - raw_fd / rawy
raa_other_old = raa_other
raa_other = (
frac_cent
* rawy
* sigma_mb
/ 2
/ acc_eff_same
/ delta_pt
/ delta_y
/ br
/ n_events
)
raa_other = frac_cent * rawy * sigma_mb / 2 / acc_eff_same / delta_pt / delta_y / br / n_events
delta_raa = abs((raa_other - raa_other_old) / raa_other)
else:
if raa_rat == 1.0 and taa == 1.0: # pp
frac.append(
1
- sigma
* delta_pt
* delta_y
* acc_eff_other
* br
* n_events
* 2
/ rawy
/ sigma_mb
)
frac.append(1 - sigma * delta_pt * delta_y * acc_eff_other * br * n_events * 2 / rawy / sigma_mb)
else: # p-Pb or Pb-Pb: iterative evaluation of Raa needed
delta_raa = 1.0
frac_tmp = 1.0
while delta_raa > 1.0e-3:
raw_fd = (
taa
* raa_rat
* raa_other
* sigma
* delta_pt
* delta_y
* acc_eff_other
* br
* n_events
* 2
taa * raa_rat * raa_other * sigma * delta_pt * delta_y * acc_eff_other * br * n_events * 2
)
frac_tmp = 1 - raw_fd / rawy
raa_other_old = raa_other
raa_other = (
frac_tmp
* rawy
* sigma_mb
/ 2
/ acc_eff_same
/ delta_pt
/ delta_y
/ br
/ n_events
)
raa_other = frac_tmp * rawy * sigma_mb / 2 / acc_eff_same / delta_pt / delta_y / br / n_events
delta_raa = abs((raa_other - raa_other_old) / raa_other)
frac.append(frac_tmp)

Expand Down Expand Up @@ -323,8 +233,6 @@ def get_hist_binlimits(histo):
n_limits = histo.GetNbinsX() + 1
low_edge = histo.GetBinLowEdge(1)
bin_width = histo.GetBinWidth(1)
bin_limits = np.array(
[low_edge + i_bin * bin_width for i_bin in range(n_limits)], "d"
)
bin_limits = np.array([low_edge + i_bin * bin_width for i_bin in range(n_limits)], "d")

return bin_limits
Loading

0 comments on commit e83093f

Please sign in to comment.