Skip to content
Snippets Groups Projects
Commit 03114fb7 authored by Markus Prim's avatar Markus Prim
Browse files

Add Hybrid class with methods to calculate the hybrid weights from a given MC...

Add Hybrid class with methods to calculate the hybrid weights from a given MC sample and to get the weight for a given phase space point.
parent 9b18c235
Branches
No related tags found
No related merge requests found
{
"BINS_MX": [0.0, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5],
"BINS_ELB":[0.0, 0.5, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 3.0],
"BINS_Q2": [0.0, 2.5, 5.0, 7.5, 10.0, 12.5, 15.0, 20.0, 25.0]
}
\ No newline at end of file
import numpy
import pandas
import json
import os
from pathlib import Path
BINS_MX = numpy.array([0., 1.4, 1.6, 1.8, 2., 2.5, 3., 3.5])
BINS_ELB = numpy.array([0., 0.5, 1., 1.25, 1.5, 1.75, 2., 2.25, 3.])
BINS_Q2 = numpy.array([0., 2.5, 5., 7.5, 10., 12.5, 15., 20., 25.])
class Hybrid:
def __init__(self, hybrid_config: str = None) -> None:
"""Initialize the hybrid weight with the given configuration.
The default binning is given by
* BINS_MX: [0.0, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5]
* BINS_EL":[0.0, 0.5, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 3.0]
* BINS_Q": [0.0, 2.5, 5.0, 7.5, 10.0, 12.5, 15.0, 20.0, 25.0]
which is motivated by theory. This should only be changed with proper theoretical motivation.
Parameters
----------
hybrid_config: str
Path to json file defining the binning of the Hybrid model. Has to contain the three attributes:
* BINS_MX: The binning for the Hybrid model in mX.
* BINS_ELB: The binning of the Hybrid model in El_B.
* BINS_Q2: The binning of the Hybrid model in q2.
"""
if hybrid_config is None:
this_dir = Path(__file__).parent
hybrid_config = this_dir / 'hybrid_binning.json'
with open(hybrid_config) as json_file:
config = json.load(json_file)
self.bins_mX = config['BINS_MX']
self.bins_El_B = config['BINS_ELB']
self.bins_q2 = config['BINS_Q2']
self.range_mX = (self.bins_mX[0], self.bins_mX[1])
self.range_El_B = (self.bins_El_B[0], self.bins_El_B[1])
self.range_q2 = (self.bins_q2[0], self.bins_q2[1])
def calculate_weight(self, x: numpy.array, hybrid_weights: numpy.array) -> float:
"""Calculate the weight in the 3D phase space of the inclusive model.
Parameters
----------
x : numpy.array
Phase space location of the inclusive prediction in the form (El_B, q2, mX).
hybrid_weights: numpy:array
Hybrid weights generated by :func:`~hybrid.Hybrid.generate_hybrid_weights`.
Returns
-------
hybrid_weight: float
The hybrid weight w_i at the given phase space point x.
"""
try:
return hybrid_weights[
int(numpy.digitize(x[0], self.bins_El_B)) - 1,
int(numpy.digitize(x[1], self.bins_q2)) - 1,
int(numpy.digitize(x[2], self.bins_mX)) - 1
]
except IndexError:
return 0
def generate_hybrid_weights(self, inclusive: pandas.DataFrame, exclusive: pandas.DataFrame) -> numpy.array:
"""Calculate the Hybrid weights w_i, so that the bin content in the Hybrid model is given by
H_i = R_i + w_i I_i, where H_i is the bin content of the inclusive prediction in the Hybrid model, R_i is the bin
content of the resonant contributions, and I_i is the bin content of the total inclusive prediction.
The required columns in the inclusive and exclusive data frames are:
* mX: the invariant mass of the hadronic system (the mass of the resonance).
* El_B: the lepton momentum in the B reference frame.
* q2: the momentum transfer to the lepton-neutrino system.
Parameters
----------
inclusive: pandas.DataFrame
Inclusive MC.
exclusive: pandas.DataFrame
Exclusive MC. Can contain any number of resonances.
Returns
-------
hybrid_weights: numpy.array
3D histogram containing the weights of the Hybrid model. To be used with
:func:`~hybrid.Hybrid.calculate_weight`.
"""
H_exc, _ = numpy.histogramdd(
exclusive[['El_B', 'q2', 'mX']].values,
bins=[self.bins_El_B, self.bins_q2, self.bins_mX],
range=[self.range_El_B, self.range_q2, self.range_mX],
normed=False,
weights=exclusive['__weight__']
)
H_incl, _ = numpy.histogramdd(
inclusive[['El_B', 'q2', 'mX']].values,
bins=[self.bins_El_B, self.bins_q2, self.bins_mX],
range=[self.range_El_B, self.range_q2, self.range_mX],
normed=False,
weights=inclusive['__weight__']
)
hybrid_weights = (H_incl - H_exc) / H_incl
# Replace nan values with 1
hybrid_weights = numpy.where(
numpy.isnan(hybrid_weights), numpy.ones(hybrid_weights.shape), hybrid_weights)
# Replace < 0 values with 0
hybrid_weights = numpy.where(
hybrid_weights < 0, numpy.zeros(hybrid_weights.shape), hybrid_weights)
return hybrid_weights
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment