samples.h File Reference

#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/numpy.h>
#include "Samples/SampleHandlerBase.h"
#include "TH1.h"
#include "TH2.h"

Include dependency graph for samples.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
class	PySampleHandlerInterface
	EW: As SampleHandlerBase is an abstract base class we have to do some gymnastics to get it to get it into python. More...
class	PySampleHandlerBase
	As SampleHandlerBase is an abstract base class we have to do some gymnastics to get it to get it into python. More...

Functions
std::tuple< py::array_t< M3::float_t >, py::array_t< M3::float_t > >	TH1ToNumpy (const TH1 *hist)
std::tuple< py::array_t< M3::float_t >, py::array_t< M3::float_t >, py::array_t< M3::float_t > >	TH2ToNumpy (const TH2 *hist)
void	initSamplesModule (py::module &m_samples)

Detailed Description

Author

Ewan Miller

Definition in file samples.h.

Function Documentation

initSamplesModule()

void initSamplesModule

(

py::module &

m_samples

)

Todo:

Deal with non uniform binning

Todo:

Deal with higher dimensions MaCh3 returns flattened bins, will need to figure out how to pass bin edge info into python

Todo:

Deal with non uniform binning

Todo:

Deal with non uniform binning

Todo:

Deal with higher dimensions MaCh3 returns flattened bins, will need to figure out how to pass bin edge info into python

Todo:

Deal with non uniform binning

Todo:

Deal with higher dimensions MaCh3 returns flattened bins, will need to figure out how to pass bin edge info into python

Definition at line 423 of file samples.h.

                                           {
 
    // Bind the systematic type enum that lets us set different types of systematics
    py::enum_<TestStatistic>(m_samples, "TestStatistic")
        .value("Poisson", TestStatistic::kPoisson)
        .value("Barlow_Beeston", TestStatistic::kBarlowBeeston)
        .value("Ice_Cube", TestStatistic::kIceCube)
        .value("Pearson", TestStatistic::kPearson)
        .value("Dembinski_Abdelmottele", TestStatistic::kDembinskiAbdelmotteleb)
        .value("N_Test_Statistics", TestStatistic::kNTestStatistics);
 
    py::class_<SampleHandlerInterface, PySampleHandlerInterface /* <--- trampoline*/>(m_samples, "SampleHandlerInterface")
        .def(py::init())
        
        .def(
            "reweight", 
            &SampleHandlerInterface::Reweight,
            "reweight the MC events in this sample. You will need to override this."
        )
        
        .def(
            "get_likelihood", 
            &SampleHandlerInterface::GetLikelihood,
            "Get the sample likelihood at the current point in your model space. You will need to override this."
        )
        
        .def(
            "set_test_stat",
            &SampleHandlerInterface::SetTestStatistic,
            "Set the test statistic that should be used when calculating likelihoods. \n\
            :param test_stat: The new test statistic to use",
            py::arg("test_stat")
        )
 
        .def(
            "get_bin_LLH",
            py::overload_cast<double, double, double>(&SampleHandlerInterface::GetTestStatLLH, py::const_),
            "Get the LLH for a bin by comparing the data and MC. The result depends on having previously set the test statistic using :py:meth:`pyMaCh3.samples.SampleHandlerInterface.set_test_stat` \n\
            :param data: The data content of the bin. \n\
            :param mc: The mc content of the bin \n\
            :param w2: The Sum(w_{i}^2) (sum of weights squared) in the bin, which is sigma^2_{MC stats}",
            py::arg("data"), 
            py::arg("mc"), 
            py::arg("w2")
        )
    ; // End of SampleHandlerInterface binding
 
    py::class_<SampleHandlerBase, PySampleHandlerBase /* <--- trampoline*/, 
            SampleHandlerInterface>(m_samples, "SampleHandlerBase")
        .def(
            py::init<std::string, ParameterHandlerGeneric*>(),
            "This should never be called directly as SampleHandlerBase is an abstract base class. \n\
            However when creating a derived class, in the __init__() method, you should call the parent constructor i.e. this one by doing:: \n\
            \n\
            \tsuper(<your derived SampleHandler class>, self).__init__(*args) \n\
            \n ",
            py::arg("mc_version"), py::arg("xsec_cov"))
 
        .def(
            "get_mc_hist",
            [](SampleHandlerBase &self, const int sample) {
 
              int Dimension = self.GetNDim(sample);
 
              //self.Reweight();
              
              // Get the histogram pointer BEFORE cloning
              const TH1 *hist_original = self.GetMCHist(sample);
              
              // Debug: Check the original histogram
              if (!hist_original) {
                throw std::runtime_error("GetMCHist returned null pointer");
              }
              
              // Now clone it
              TH1D *hist = static_cast<TH1D*>(hist_original->Clone("cloned_hist"));
 
              if (Dimension == 1) {
                // 1D histogram
                auto [contents, edgesX] = TH1ToNumpy(hist);
                auto edgesY = py::array_t<M3::float_t>();
                return py::make_tuple(contents, edgesX, edgesY);
              } else if (Dimension == 2) {
 
                TH2Poly *hist2poly = dynamic_cast<TH2Poly *>(hist);
                if (hist2poly) {
                    throw std::runtime_error("pyMaCh3 can't do non-uniform binning for now :(");
                }
 
                // 2D histogram - cast to TH2
                TH2 *hist2d = dynamic_cast<TH2 *>(hist);
                if (!hist2d) {
                  throw std::runtime_error("Failed to cast to TH2");
                }
                auto [contents, edgesX, edgesY] = TH2ToNumpy(hist2d);
                return py::make_tuple(contents, edgesX, edgesY);
              } else {
                throw std::invalid_argument("Dimension must be 1 or 2");
              }
            },
            py::return_value_policy::reference_internal,
            py::arg("sample"),
            "Get MC histogram as numpy arrays.\n"
            "For 1D: Returns (contents, edges)\n"
            "For 2D: Returns (contents, edgesX, edgesY)\n"
            "where contents is shape (nbinsY, nbinsX) for 2D")
 
        .def(
            "get_data_hist",
            [](SampleHandlerBase &self, const int sample) {
              const TH1 *hist = self.GetDataHist(sample);
 
              int Dimension = self.GetNDim(sample);
              
              if (Dimension == 1) {
                // 1D histogram
                const auto [contents, edgesX] = TH1ToNumpy(hist);
                const auto edgesY = py::array_t<M3::float_t>();
                return py::make_tuple(contents, edgesX, edgesY);
              } else if (Dimension == 2) {
 
                const TH2Poly *hist2poly = dynamic_cast<const TH2Poly *>(hist);
                if (hist2poly) {
                    throw std::runtime_error("pyMaCh3 can't do non-uniform binning for now :(");
                }
                
                throw std::runtime_error("pyMaCh3 can't do non-uniform binning for now :(");
 
                // 2D histogram - cast to TH2
                const TH2 *hist2d = dynamic_cast<const TH2 *>(hist);
                if (!hist2d) {
                  throw std::runtime_error("Failed to cast to TH2");
                }
                const auto [contents, edgesX, edgesY] = TH2ToNumpy(hist2d);
                return py::make_tuple(contents, edgesX, edgesY);
              } else {
                throw std::invalid_argument("Dimension must be 1 or 2");
              }
            },
            py::arg("Dimension"),
            "Get Data histogram as numpy arrays.\n"
            "For 1D: Returns (contents, edges)\n"
            "For 2D: Returns (contents, edgesX, edgesY)\n"
            "where contents is shape (nbinsY, nbinsX) for 2D")
 
        .def(
            "get_w2_hist",
            [](SampleHandlerBase &self, const int sample) {
              const TH1 *hist = self.GetW2Hist(sample);
 
              int Dimension = self.GetNDim(sample);
 
              if (Dimension == 1) {
                // 1D histogram
                const auto [contents, edgesX] = TH1ToNumpy(hist);
                const auto edgesY = py::array_t<M3::float_t>();
                return py::make_tuple(contents, edgesX, edgesY);
              } else if (Dimension == 2) {
 
                const TH2Poly *hist2poly = dynamic_cast<const TH2Poly *>(hist);
                if (hist2poly) {
                    throw std::runtime_error("pyMaCh3 can't do non-uniform binning for now :(");
                }
                
                // 2D histogram - cast to TH2
                const TH2 *hist2d = dynamic_cast<const TH2 *>(hist);
                if (!hist2d) {
                  throw std::runtime_error("Failed to cast to TH2");
                }
                const auto [contents, edgesX, edgesY] = TH2ToNumpy(hist2d);
                return py::make_tuple(contents, edgesX, edgesY);
              } else {
                throw std::invalid_argument("Dimension must be 1 or 2");
              }
            },
            py::arg("sample"),
            "Get W2 histogram as numpy arrays.\n"
            "For 1D: Returns (contents, edges)\n"
            "For 2D: Returns (contents, edgesX, edgesY)\n"
            "where contents is shape (nbinsY, nbinsX) for 2D");
 
    /* Not sure if this will be needed in future versions of MaCh3 so leaving commented for now
    py::class_<fdmc_base>(m_samples, "MCstruct")
        .def(py::init())
        
        // Because a lot of the variables in fdmc_base use c style arrays,
        // we need to provide some setter functions to be able to set them using more
        // "pythony" objects, e.g. lists and numpy arrays
        .def(
            "set_event_variable_values", 
            [](fdmc_base &self, int dim, py::array_t<M3::float_t, py::array::c_style> &array)
            {
                py::buffer_info bufInfo = array.request();
 
                if ( dim > 2 )
                    throw MaCh3Exception(__FILE__, __LINE__, "Currently only dimensions of 1 or 2 are supported sorry :(");
 
                if ( bufInfo.ndim != 1 )
                    throw MaCh3Exception(__FILE__, __LINE__, "Number of dimensions in parameter array must be one if setting only one of the event variable arrays!");
                
                if( dim ==1 )
                    self.x_var = array.data();
                    
                else if ( dim == 2)
                    self.y_var = array.data();
            }
        )
    ;
    */
}

TH1ToNumpy()

std::tuple<py::array_t<M3::float_t>, py::array_t<M3::float_t> > TH1ToNumpy

(

const TH1 *

hist

)

Definition at line 17 of file samples.h.

                                                                                 {
    if (!hist) {
        throw std::runtime_error("Histogram pointer is null");
    }
    
    int nbins = hist->GetNbinsX();
    
    // Create numpy array for bin contents
    py::array_t<M3::float_t> contents(nbins);
    auto contents_buf = contents.request();
    M3::float_t* contents_ptr = static_cast<M3::float_t*>(contents_buf.ptr);
    
    // Create numpy array for bin edges (nbins + 1 edges)
    py::array_t<M3::float_t> edges(nbins + 1);
    auto edges_buf = edges.request();
    M3::float_t* edges_ptr = static_cast<M3::float_t*>(edges_buf.ptr);
    
    // Copy bin contents (ROOT bins start at 1, not 0)
    for (int i = 0; i < nbins; ++i) {
        contents_ptr[i] = hist->GetBinContent(i + 1);
    }
    
    // Copy bin edges
    for (int i = 0; i <= nbins; ++i) {
        edges_ptr[i] = hist->GetBinLowEdge(i + 1);
    }
    // Add the upper edge of the last bin
    edges_ptr[nbins] = hist->GetBinLowEdge(nbins + 1) + hist->GetBinWidth(nbins + 1);
    
    return std::make_tuple(contents, edges);
}

TH2ToNumpy()

std::tuple<py::array_t<M3::float_t>, py::array_t<M3::float_t>, py::array_t<M3::float_t> > TH2ToNumpy

(

const TH2 *

hist

)

Definition at line 50 of file samples.h.

                                                                                                       {
    if (!hist) {
        throw std::runtime_error("Histogram pointer is null");
    }
    
    int nbinsX = hist->GetNbinsX();
    int nbinsY = hist->GetNbinsY();
    
    // Create 2D numpy array for bin contents (shape: nbinsY x nbinsX to match numpy convention)
    py::array_t<M3::float_t> contents({nbinsY, nbinsX});
    auto contents_buf = contents.request();
    M3::float_t* contents_ptr = static_cast<M3::float_t*>(contents_buf.ptr);
    
    // Create numpy arrays for bin edges
    py::array_t<M3::float_t> edgesX(nbinsX + 1);
    auto edgesX_buf = edgesX.request();
    M3::float_t* edgesX_ptr = static_cast<M3::float_t*>(edgesX_buf.ptr);
    
    py::array_t<M3::float_t> edgesY(nbinsY + 1);
    auto edgesY_buf = edgesY.request();
    M3::float_t* edgesY_ptr = static_cast<M3::float_t*>(edgesY_buf.ptr);
    
    // Copy bin contents (ROOT bins start at 1, not 0)
    // Note: numpy uses row-major order (C-style), so we iterate Y then X
    for (int iy = 0; iy < nbinsY; ++iy) {
        for (int ix = 0; ix < nbinsX; ++ix) {
            contents_ptr[iy * nbinsX + ix] = hist->GetBinContent(ix + 1, iy + 1);
        }
    }
    
    // Copy X bin edges
    for (int i = 0; i <= nbinsX; ++i) {
        edgesX_ptr[i] = hist->GetXaxis()->GetBinLowEdge(i + 1);
    }
    edgesX_ptr[nbinsX] = hist->GetXaxis()->GetBinLowEdge(nbinsX + 1) + 
                         hist->GetXaxis()->GetBinWidth(nbinsX + 1);
    
    // Copy Y bin edges
    for (int i = 0; i <= nbinsY; ++i) {
        edgesY_ptr[i] = hist->GetYaxis()->GetBinLowEdge(i + 1);
    }
    edgesY_ptr[nbinsY] = hist->GetYaxis()->GetBinLowEdge(nbinsY + 1) + 
                         hist->GetYaxis()->GetBinWidth(nbinsY + 1);
    
    return std::make_tuple(contents, edgesX, edgesY);
}