Skip to content
Snippets Groups Projects
Commit aa53bbcd authored by Johannes Schumann's avatar Johannes Schumann
Browse files

Resolve "Add helicity probability"

parent 2c78f4bf
No related branches found
No related tags found
1 merge request!79Resolve "Add helicity probability"
Unreleased changes
------------------
* Add field for right handed helicity probability for outgoing lepton in GiBUUOutput
v1.6.1
----------------------------
......
......@@ -464,6 +464,30 @@ class GiBUUOutput:
y = pq / pk
return y
@staticmethod
def right_helicity_probability(roottuple_data):
"""
Calculate the probability for the outgoing lepton from the neutrino vertex
to have a right handed helicity
Definition: p = ( 1-|p|/(E+m) )/2
Parameters
----------
roottuple_data: awkward.highlevel.Array
"""
sec_lepton_pdgid = np.sign(
roottuple_data.process_ID) * (7 + 2 * roottuple_data.flavor_ID +
np.abs(roottuple_data.process_ID))
sec_lepton_mass = np.array(
[Particle.from_pdgid(p).mass * 1e-3 for p in sec_lepton_pdgid])
total_lepton_P = np.sqrt(roottuple_data.lepOut_Px**2 +
roottuple_data.lepOut_Py**2 +
roottuple_data.lepOut_Pz**2)
rh_probability = 0.5 * (1 - total_lepton_P /
(roottuple_data.lepOut_E + sec_lepton_mass))
return rh_probability
@property
def A(self):
grp = self.jobcard["target"]
......@@ -539,6 +563,7 @@ class GiBUUOutput:
visEfrac = visible_energy_fraction(ak.flatten(retval.E),
ak.flatten(retval.barcode))
retval["visEfrac"] = ak.unflatten(visEfrac, counts)
retval["rh_prob"] = GiBUUOutput.right_helicity_probability(retval)
return retval
@property
......
......@@ -106,6 +106,11 @@ class TestGiBUUOutput(unittest.TestCase):
mask = self.output.free_particle_mask
np.testing.assert_array_equal(mask[0], [True, False, True, True, True])
def test_right_handed_helicity(self):
arr = self.output.arrays
np.testing.assert_array_almost_equal(
arr.rh_prob[:3], [4.518065e-05, 9.428467e-06, 3.488926e-05])
@pytest.mark.skipif(not KM3NET_LIB_AVAILABLE,
reason="KM3NeT dataformat required")
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment