# ===================================================
#
# Copyright (c) 2023-2024
# SPARKX Team
#
# GNU General Public License (GPLv3 or later)
#
# ===================================================
import numpy as np
import warnings
from typing import Union, List
[docs]
class CentralityClasses:
"""
Class for defining centrality classes based on event multiplicity.
.. note::
It is the user's responsibility to ensure that the amount of events used
to determine the centrality classes is sufficient to provide reliable
results. The recommended minimum number of events is at least a few
hundred.
Parameters
----------
events_multiplicity : list or numpy.ndarray
List or array containing the multiplicity values for each event.
centrality_bins : list or numpy.ndarray
List or array defining the boundaries of centrality classes as percentages.
Raises
------
TypeError
If :code:`events_multiplicity` or :code:`centrality_bins` is not a list
or :code:`numpy.ndarray`.
Attributes
----------
events_multiplicity_ : list or numpy.ndarray
Stores the input multiplicity values for each event.
centrality_bins_ : list or numpy.ndarray
Stores the input boundaries of centrality classes.
dNchdetaMin_ : list
Minimum values of multiplicity for each centrality class.
dNchdetaMax_ : list
Maximum values of multiplicity for each centrality class.
dNchdetaAvg_ : list
Average values of multiplicity for each centrality class.
dNchdetaAvgErr_ : list
Average errors of multiplicity for each centrality class.
Methods
-------
get_centrality_class:
Return the index of the centrality bin for a given multiplicity value.
output_centrality_classes:
Write centrality class information to a file.
Examples
--------
.. highlight:: python
.. code-block:: python
:linenos:
>>> from sparkx import *
>>> import random as rd
>>> rd.seed(0)
>>> # generate 300 random event multiplicities between 50 and 1500
>>> events_multiplicity = [rd.randint(50, 1500) for _ in range(300)]
>>> centrality_bins = [0, 10, 30, 50, 70, 90, 100]
>>> centrality_obj = CentralityClasses(events_multiplicity=events_multiplicity,
... centrality_bins=centrality_bins)
>>> centrality_obj.get_centrality_class(1490)
0
>>> centrality_obj.output_centrality_classes('centrality_output.txt')
"""
def __init__(
self,
events_multiplicity: Union[List[float], np.ndarray],
centrality_bins: Union[List[float], np.ndarray],
) -> None:
if not isinstance(events_multiplicity, (list, np.ndarray)):
raise TypeError(
"'events_multiplicity' is not list or numpy.ndarray"
)
if not isinstance(centrality_bins, (list, np.ndarray)):
raise TypeError("'centrality_bins' is not list or numpy.ndarray")
# Check if centrality_bins is sorted
if not all(
centrality_bins[i] <= centrality_bins[i + 1]
for i in range(len(centrality_bins) - 1)
):
warnings.warn(
"'centrality_bins' is not sorted. Sorting automatically."
)
centrality_bins.sort()
# Check for uniqueness of values
# Remove duplicates from the list
unique_bins = []
seen = set()
multiple_same_entries = False
for item in centrality_bins:
if item not in seen:
unique_bins.append(item)
seen.add(item)
else:
multiple_same_entries = True
if multiple_same_entries:
warnings.warn(
"'centrality_bins' contains duplicate values. They are removed automatically."
)
# Check for negative values and values greater than 100
if any(value < 0.0 or value > 100.0 for value in centrality_bins):
raise ValueError(
"'centrality_bins' contains values less than 0 or greater than 100."
)
self.events_multiplicity_ = events_multiplicity
self.centrality_bins_ = unique_bins
self.dNchdetaMin_: List[float] = []
self.dNchdetaMax_: List[float] = []
self.dNchdetaAvg_: List[float] = []
self.dNchdetaAvgErr_: List[float] = []
self.__create_centrality_classes()
def __create_centrality_classes(self) -> None:
"""
Create centrality classes based on event multiplicity.
Parameters
----------
None
Returns
-------
None
Raises
------
ValueError
If the number of events is less than 4.
If the multiplicity in 'events_multiplicity' is negative.
Notes
-----
This function creates four sub-samples for error determination based on the
event multiplicity. The numbers are distributed evenly to the four sub-samples,
which are then sorted in descending order. The averages and errors are
calculated for each sub-sample. Finally, the events are sorted by multiplicity,
and boundaries of centrality classes are determined.
"""
number_events = len(self.events_multiplicity_)
if number_events < 4:
raise ValueError("The number of events has to be larger than 4")
# create four sub samples for the error determination
event_sample_A = []
event_sample_B = []
event_sample_C = []
event_sample_D = []
# distribute the numbers evenly
for i, multiplicity in enumerate(self.events_multiplicity_):
if multiplicity < 0:
raise ValueError(
"Multiplicity in 'events_multiplicity' is negative"
)
if i % 4 == 0:
event_sample_A.append(multiplicity)
elif i % 4 == 1:
event_sample_B.append(multiplicity)
elif i % 4 == 2:
event_sample_C.append(multiplicity)
elif i % 4 == 3:
event_sample_D.append(multiplicity)
event_sample_A = sorted(event_sample_A, reverse=True)
event_sample_B = sorted(event_sample_B, reverse=True)
event_sample_C = sorted(event_sample_C, reverse=True)
event_sample_D = sorted(event_sample_D, reverse=True)
MinRecord = int(number_events / 4 * self.centrality_bins_[0] / 100.0)
for i in range(1, len(self.centrality_bins_)):
MaxRecord = int(
number_events / 4 * self.centrality_bins_[i] / 100.0
)
AvgA = np.mean(event_sample_A[MinRecord:MaxRecord])
AvgB = np.mean(event_sample_B[MinRecord:MaxRecord])
AvgC = np.mean(event_sample_C[MinRecord:MaxRecord])
AvgD = np.mean(event_sample_D[MinRecord:MaxRecord])
Avg = (AvgA + AvgB + AvgC + AvgD) / 4.0
Err = np.sqrt(
(
(AvgA - Avg) ** 2
+ (AvgB - Avg) ** 2
+ (AvgC - Avg) ** 2
+ (AvgD - Avg) ** 2
)
/ 3.0
)
self.dNchdetaAvg_.append(Avg)
self.dNchdetaAvgErr_.append(Err)
MinRecord = MaxRecord
# sort events by multiplicity and determine boundaries of centrality
# classes
global_event_record = sorted(self.events_multiplicity_, reverse=True)
MinRecord = int(number_events * self.centrality_bins_[0] / 100.0)
for i in range(1, len(self.centrality_bins_)):
MaxRecord = int(number_events * self.centrality_bins_[i] / 100.0)
self.dNchdetaMax_.append(global_event_record[MinRecord])
self.dNchdetaMin_.append(global_event_record[MaxRecord - 1])
MinRecord = MaxRecord
[docs]
def get_centrality_class(self, dNchdEta: float) -> int:
"""
This function determines the index of the centrality bin for a given
multiplicity value based on the predefined centrality classes.
In the case that the multiplicity input exceeds the largest or smallest
value of the multiplicity used to determine the centrality classes, the
function returns the index of the most central or most peripheral bin,
respectively.
Parameters
----------
dNchdEta : float
Multiplicity value.
Returns
-------
int
Index of the centrality bin.
"""
# check if the multiplicity is in the most central bin
if dNchdEta >= self.dNchdetaMin_[0]:
return 0
# check if the multiplicity is in the most peripheral bin
elif dNchdEta < self.dNchdetaMin_[len(self.dNchdetaMin_) - 2]:
return len(self.dNchdetaMin_) - 1
# check if the multiplicity is in one of the intermediate bins
else:
for i in range(1, len(self.dNchdetaMin_) - 1):
if (dNchdEta >= self.dNchdetaMin_[i]) and (
dNchdEta < self.dNchdetaMin_[i - 1]
):
return i
# default case to satisfy mypy's requirement for a return statement
return -1
[docs]
def output_centrality_classes(self, fname: str) -> None:
"""
Write centrality class information to a file.
Parameters
----------
fname : str
Name of the output file.
Raises
------
TypeError
If :code:`fname` is not a string.
Notes
-----
This function writes the centrality class information, including minimum,
maximum, average multiplicities, and average errors, to the specified file.
"""
# Check if fname is a string
if not isinstance(fname, str):
raise TypeError("'fname' should be a string.")
# Write the information to the file
with open(fname, "w") as out_stream:
out_stream.write(
"# CentralityMin CentralityMax dNchdEtaMin dNchdEtaMax dNchdEtaAvg dNchdEtaAvgErr\n"
)
for i in range(1, len(self.dNchdetaMin_)):
out_stream.write(
f"{self.centrality_bins_[i - 1]} - {self.centrality_bins_[i]} "
f"{self.dNchdetaMin_[i - 1]} {self.dNchdetaMax_[i - 1]} "
f"{self.dNchdetaAvg_[i - 1]} {self.dNchdetaAvgErr_[i - 1]}\n"
)