Merge pull request #6 from MadAnalysis/histos

Histogram Reader

README.md

@@ -8,7 +8,12 @@
 ## Installation
 `pip install ma5-expert`
 
-## Cutflow Collection
+## Outline
+* [Cutflow Collection](#cutflow-collection)
+* [Histogram Collection](#histogram-collection)
+* [Citation](#citation)
+
+### Cutflow Collection
 
  * Parse all the signal regions and construct an object-base, interactable cutflow.
  * Write combined LaTeX tables for different samples.
@@ -141,7 +146,41 @@ ATLAS.addSignalRegion('SRA_H', ma5.SRA_H.CutNames, SRA_presel+[7.0])
 ```
 where all properties shown above applies to this new object as well.
 
-## Citation 
+### Histogram Collection
+
+* Parse all the histograms available in the `Histos.saf` file into interactable histogram object.
+
+```python
+import ma5_expert as ma5
+collection = ma5.histogram.Collection(
+    "examples/mass1000005_300.0_mass1000022_60.0_mass1000023_250.0_xs_5.689/Output/SAF/defaultset/atlas_susy_2018_31/Histograms/histos.saf", 
+    xsection=5.689, lumi=139.
+)
+
+print(collection)
+# Collection of 6 histograms from examples/mass1000005_300.0_mass1000022_60.0_mass1000023_250.0_xs_5.689/Output/SAF/defaultset/atlas_susy_2018_31/Histograms
+#    * MadAnalysis 5 Histogram: SRA_Meff
+#    * MadAnalysis 5 Histogram: SRA_Mh
+#    * MadAnalysis 5 Histogram: SRB_PTj1
+#    * MadAnalysis 5 Histogram: SRB_MhAvg
+#    * MadAnalysis 5 Histogram: SRC_MET
+#    * MadAnalysis 5 Histogram: SRC_Sig
+```
+
+Extract the plotting information:
+```python
+xbins, bins, weights = collection.lumi_histogram("SRA_Mh")
+plt.hist(xbins, bins=bins, weights=weights)
+plt.xlabel("$M_{h}\ {\\rm [GeV]}$")
+plt.ylabel("${\\rm Number\ of\ events}$")
+plt.xlim([min(bins), max(bins)])
+plt.show()
+```
+<p align="center">
+<img src="./examples/SRA_Mh.png" alt="SRA_Mh" style="width:400px;"/>
+</p>
+
+### Citation 
 Developed for [arXiv:2006.09387](http://arxiv.org/abs/2006.09387)
 
 ```bibtex

setup.py

@@ -12,15 +12,15 @@
 
 setup(
     name="ma5_expert",
-    version="1.0.5",
+    version="2.0.0",
     description=("MadAnalysis 5 interpreter for Expert mode"),
     long_description=long_description,
     long_description_content_type="text/markdown",
     url="https://github.com/MadAnalysis/ma5_expert",
     project_urls={
         "Bug Tracker": "https://github.com/MadAnalysis/ma5_expert/issues",
     },
-    download_url = "https://github.com/MadAnalysis/ma5_expert/archive/refs/tags/v1.0.4.tar.gz",
+    download_url = "https://github.com/MadAnalysis/ma5_expert/archive/refs/tags/v2.0.0.tar.gz",
     author="Jack Y. Araz",
     author_email=("[email protected]"),
     license="MIT",

src/ma5_expert/__init__.py

@@ -10,5 +10,6 @@
 log.setLevel(logging.INFO)
 
 from ma5_expert import cutflow
+from ma5_expert import histogram
 
 __all__ = cutflow.__all__

src/ma5_expert/histogram/__init__.py

@@ -0,0 +1 @@
+from .reader import Collection

src/ma5_expert/histogram/bin.py

@@ -0,0 +1,45 @@
+class Bin:
+    """
+    Histogram bin object
+
+    Parameters
+    ----------
+    sumw: float
+        sum of weights
+    bin_max: float
+        maximum bin limit
+    bin_min: float
+        minimum bin limit
+    underflow: bool
+        is this underflow bin
+    overflow: bool
+        is this overflow bin
+    """
+    def __init__(
+        self,
+        sumw: float,
+        underflow: bool = False,
+        overflow: bool = False,
+    ):
+        self._sumw = sumw
+        self._underflow = underflow
+        self._overflow = overflow
+
+    @property
+    def isUnderflow(self) -> bool:
+        """ Is this underflow bin? """
+        return self._underflow
+
+    @property
+    def isOverflow(self) -> bool:
+        """ Is this overflow bin? """
+        return self._overflow
+
+    @property
+    def sumW(self) -> float:
+        """ Sum of Weights """
+        return self._sumw
+
+    def eff(self, total_sumw: float) -> float:
+        """ Return bin efficiency """
+        return self._sumw / total_sumw if total_sumw > 0 else float("inf")

src/ma5_expert/histogram/histo.py

@@ -0,0 +1,156 @@
+import numpy as np
+
+from .bin import Bin
+from typing import Sequence, Text, Union, MutableSequence
+
+
+class Histo:
+    """
+    Object oriented Histogram definition
+    """
+
+    name: Text
+    ID: int
+    _nbins: int
+    regions: Sequence[Text]
+    _nEvents: int
+    _normEwEvents: float
+    _nEntries: float
+    _normEwEntries: float
+    _sumWeightsSq: float
+    _sumValWeight: float
+    _sumValSqWeight: float
+    _xmin: float
+    _xmax: float
+    _bins: MutableSequence[Bin]
+    _normalisation_frac: Union[float, Text]
+
+    def __init__(self):
+        self.ID = -1
+        self.name = "__unknown_histo__"
+        self._nbins = 0
+        self.regions = []
+        self._nEvents = 0
+        self._normEwEvents = 0
+        self._nEntries = 0
+        self._normEwEntries = 0
+        self._sumWeightsSq = 0
+        self._sumValWeight = 0
+        self._sumValSqWeight = 0
+        self._xmin = 0
+        self._xmax = 0
+        self._bins = []
+        self._normalisation_frac = "_normEwEvents"
+
+    def __repr__(self):
+        return f"MadAnalysis 5 Histogram: {self.name}"
+
+    @property
+    def weight_normalisation(self) -> float:
+        """
+        Retrun the weight normalisation factor. Typically this is the value to normalise the
+        sum of weights per bin, hence the total sum of weights of the sample. By default,
+        weigted sum of entries are used.
+        """
+        if isinstance(self._normalisation_frac, str):
+            return getattr(self, self._normalisation_frac)
+        else:
+            return self._normalisation_frac
+
+    @weight_normalisation.setter
+    def weight_normalisation(self, value: float):
+        if not isinstance(value, float):
+            raise ValueError("Input can only be float.")
+        self._normalisation_frac = value
+
+    @property
+    def size(self):
+        return self._nbins
+
+    def _add_bin(self, bin_info: dict) -> None:
+        """
+        Add a bin to the histogram
+
+        Parameters
+        ----------
+        bin_info: dict
+            bin information
+
+        Raises
+        -------
+        AssertionError:
+            If the histogram name and ID information does not match the info can not be added.
+        """
+        if self.ID == -1 and self.name == "__unknown_histo__":
+            self.ID = int(bin_info["ID"])
+            self.name = bin_info["name"]
+            self._nbins = int(bin_info["nbins"])
+            self.regions = bin_info["region"]
+            self._nEvents = int(bin_info["nEvents"])
+            self._normEwEvents = float(bin_info["normEwEvents"])
+            self._nEntries = float(bin_info["normEwEvents"])
+            self._normEwEntries = float(bin_info["normEwEntries"])
+            self._sumWeightsSq = float(bin_info["sumWeightsSq"])
+            self._sumValWeight = float(bin_info["sumValWeight"])
+            self._sumValSqWeight = float(bin_info["sumValSqWeight"])
+            self._xmin = float(bin_info["xmin"])
+            self._xmax = float(bin_info["xmax"])
+            self._bins.append(
+                Bin(float(bin_info["value"]), bin_info["isUnderflow"], bin_info["isOverflow"])
+            )
+
+        else:
+            assert (
+                self.ID == int(bin_info["ID"]) and self.name == bin_info["name"]
+            ), "Merging different types of histograms are not allowed."
+            self._bins.append(
+                Bin(float(bin_info["value"]), bin_info["isUnderflow"], bin_info["isOverflow"])
+            )
+
+    def _w(self, weight: float):
+        return np.array(
+            [
+                b.eff(self.weight_normalisation) * weight
+                for b in self._bins
+                if not (b.isOverflow or b.isUnderflow)
+            ],
+            dtype=np.float32,
+        )
+
+    @property
+    def weights(self) -> np.ndarray:
+        return self._w(1.0)
+
+    def norm_weights(self, xsec: float) -> np.ndarray:
+        """
+        Normalised bin weights with respect to cross-section.
+        This function does not include overflow or underflow bins
+
+        Parameters
+        ----------
+        xsec: float
+            cross section in pb
+        """
+        return self._w(xsec)
+
+    def lumi_weights(self, xsec: float, lumi: float) -> np.ndarray:
+        """
+        Normalised bin weights with respect to cross section and luminosity.
+        This function does not include overflow or underflow bins
+
+        Parameters
+        ----------
+        xsec: float
+            cross-section in pb
+        lumi: float
+            luminosity in 1/fb
+        """
+        return self._w(xsec * 1000.0 * lumi)
+
+    @property
+    def bins(self) -> np.ndarray:
+        return np.linspace(self._xmin, self._xmax, self.size + 1)
+
+    @property
+    def xbins(self) -> np.ndarray:
+        return np.linspace(self._xmin, self._xmax, self.size)