andrewfowlie / lhco_reader Goto Github PK

View Code? Open in Web Editor NEW

3.0 3.0 4.0 1.23 MB

Read LHCO files into a Python class

Home Page: http://lhco-reader.readthedocs.org

License: GNU General Public License v2.0

Python 100.00%

lhco_reader's Introduction

Hi there 👋

I'm a physicist working at XJTLU, Suzhou, China 🇨🇳
🔭 I’m currently working on particle physics, particularly phase transitions and gravitational waves, and Bayesian statistics and computation
👯 I’m looking to collaborate on more or less anything. If you have an idea, feel free to write me at [email protected]
📰 Take a look at my papers
📣 and at my talks

lhco_reader's People

Contributors

Stargazers

Watchers

Forkers

lilymag adarshp fedelamagna anafoguel

lhco_reader's Issues

Error when printing events: TypeError: init() got an unexpected keyword argument 'header'

Hello,

print events[0] works well, but print events does not work:

    Python 2.7.3 (default, Apr 20 2012, 22:44:07) 
    [GCC 4.6.3] on linux2
    Type "help", "copyright", "credits" or "license" for more information.
    >>> import LHCO_reader
    >>> events = LHCO_reader.Events(f_name="example.lhco")
    >>> print events[0]
    +----------+--------+-------+--------+-------+------+------+-------+
    |  Object  |  eta   |  phi  |   PT   | jmass | ntrk | btag | hadem |
    +----------+--------+-------+--------+-------+------+------+-------+
    | electron | -0.745 | 4.253 | 286.72 |  0.0  | -1.0 | 0.0  |  0.0  |
    | electron | -0.073 | 4.681 | 44.56  |  0.0  | 1.0  | 0.0  |  0.0  |
    |   jet    | -0.565 | 1.126 | 157.44 | 12.54 | 16.0 | 0.0  |  0.57 |
    |   jet    | -0.19  | 1.328 | 130.96 |  12.3 | 18.0 | 0.0  | 10.67 |
    |   jet    | 0.811  | 6.028 | 17.49  |  3.47 | 8.0  | 0.0  |  2.37 |
    |   jet    | 0.596  | 0.853 | 12.47  |  2.53 | 7.0  | 0.0  |  1.26 |
    |   jet    | -1.816 | 0.032 |  6.11  |  1.18 | 0.0  | 0.0  |  0.56 |
    |   jet    | 0.508  | 1.421 |  6.01  |  0.94 | 7.0  | 0.0  |  2.59 |
    |   MET    |  0.0   | 2.695 | 21.43  |  0.0  | 0.0  | 0.0  |  0.0  |
    +----------+--------+-------+--------+-------+------+------+-------+
    >>> print events
    Traceback (most recent call last):
      File "<stdin>", line 1, in <module>
      File "LHCO_reader.py", line 463, in __str__
        table = pt(header=False)
    TypeError: __init__() got an unexpected keyword argument 'header'

Function for returning non-$b$-tagged jets

The b_jets() function returns $b$-tagged jets. There should be a function for non-$b$-tagged jets too.

I think the best way is that the b_jets() function takes an optional keyword True by default, tagged=True, which if False returns non-$b$-tagged jets rather than $b$-tagged jets:

events[0].b_tags()  # Returns $b$-tagged jets
events[0].b_tags(tagged=False)   # Returns non-$b$-tagged jets

Bug in $\Delta R$ function

Dear Andrew,

In the b98d35a, delta_R is defined here as

delta_R = ((o1["eta"] - o2["eta"])**2 + (o1["phi"] - o2["phi"])**2)**0.5

However, (o1["phi"] - o2["phi"]) has to be taken between $-\pi$ and $\pi$. The function of delta_R does not consider this limit, and will give the wrong results for (o1["phi"] - o2["phi"]) > 3.1415 or (o1["phi"] - o2["phi"]) < -3.1415.

Bug in $\gamma$-function

Hi Andrew,

The gamma-function in the LHCO_reader is

gamma = self[0] / abs(self)  # gamma = E / M
return gamma

Due to the uncertainty of computer, M² = E² - Px²- Py²- Pz² sometimes equals to 0 or smaller than 0. Then gamma will return the wrong results.

Fox example, for the two events in the attached file, the gamma-function will return the unphysical results for tau. In fact, the two events are taken from another file in which the gamma-function may return the correct results of tau for other events since M² of tau in these events is larger than 0.

from LHCO_reader import LHCO_reader

events = LHCO_reader.Events(f_name="gammatau.lhco")
ta0 = events[0]["tau"][0]
ta1 = events[1]["tau"][0]

print(ta0.vector().gamma(),abs(ta0.vector()))
# LHCO_reader/LHCO_reader.py:2122: RuntimeWarning: divide by zero encountered in double_scalars
#  gamma = self[0] / abs(self)  # gamma = E / M
# (inf, 0.0)

print(ta0.vector()[0]*ta0.vector()[0]-ta0.vector()[1]*ta0.vector()[1]-ta0.vector()[2]*ta0.vector()[2]-ta0.vector()[3]*ta0.vector()[3])
#0.0

print(ta1.vector().gamma(),abs(ta1.vector()))
# LHCO_reader/LHCO_reader.py:1957: RuntimeWarning: invalid value encountered in double_scalars
#  return self.__mul__(self)**0.5
# (nan, nan)

print(ta1.vector()[0]*ta1.vector()[0]-ta1.vector()[1]*ta1.vector()[1]-ta1.vector()[2]*ta1.vector()[2]-ta1.vector()[3]*ta1.vector()[3])
# -1.09139364213e-11

File

# LHCO file created with LHCO_reader (https://github.com/innisfree/LHCO_reader).
# See http://madgraph.phys.ucl.ac.be/Manual/lhco.html for a description of the LHCO format.

# +------------------+--------------------+
# | Number of events | 2                  |
# | Description      | Hand-picked events |
# +------------------+--------------------+

# Event number: 0
#          number     trig word 
0          9170       0         
#          typ        eta        phi        pt         jmass      ntrk       btag       had/em     dummy      dummy     
1          2          0.855      1.862      158.28     0.11       1.0        6.0        158.0      0.0        0.0       
2          2          0.445      1.008      59.83      0.11       -1.0       6.0        60.0       0.0        0.0       
3          2          -0.761     -0.43      55.63      0.11       1.0        6.0        56.0       0.0        0.0       
4          2          -0.072     2.43       33.31      0.11       -1.0       6.0        33.07      0.0        0.0       
5          3          2.11       1.951      94.02      -0.0       -1.0       0.0        999.9      0.0        0.0       
6          4          1.543      -1.295     220.6      31.65      12.0       1.0        0.03       0.0        0.0       
7          4          2.099      -1.427     93.82      9.86       13.0       1.0        0.11       0.0        0.0       
8          6          0.0        2.654      31.78      0.0        0.0        0.0        0.0        0.0        0.0       
# Event number: 1
#          number     trig word 
0          9183       0         
#          typ        eta        phi        pt         jmass      ntrk       btag       had/em     dummy      dummy     
1          2          0.12       -2.784     452.08     0.11       -1.0       7.0        452.0      0.0        0.0       
2          2          -2.322     -3.022     31.06      0.11       -1.0       7.0        31.09      0.0        0.0       
3          2          0.823      3.085      28.3       0.11       1.0        7.0        30.0       0.0        0.0       
4          2          -1.093     1.778      26.79      0.11       1.0        7.0        27.0       0.0        0.0       
6          3          -2.606     0.665      25.02      0.0        -1.0       0.0        0.0        0.0        0.0       
7          4          -1.327     -2.451     85.74      10.13      15.0       1.0        1.88       0.0        0.0       
8          4          -2.331     -1.199     47.36      11.16      4.0        1.0        0.61       0.0        0.0       
9          6          0.0        0.289      425.03     0.0        0.0        0.0        0.0        0.0        0.0

ImportError: No module named matplotlib.pyplot

Hi,
I am very interested in your LHCO_reader. However, When I run it, a error is reported.

$ python
Python 2.7.3 (default, Apr 20 2012, 22:44:07) 

[GCC 4.6.3] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> from LHCO_reader import LHCO_reader
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "LHCO_reader/LHCO_reader.py", line 84, in <module>
    import matplotlib.pyplot as plt
ImportError: No module named matplotlib.pyplot

In addition, could you share a completed code which is used to analyze a specific signal or background.
It will be very helpful for me to understand and use your great LHCO_reader.

Many thanks

Design: leptons vs anti-leptons, jets vs b-jets

Current structure

in the current structure of an Events() class, leptons, jets and photons are fundamental detector-level objects. Leptons have a property that indicates their charge (ntrk) and jets have a property (btag) that indicates whether they were $b$-tagged.

To deal with cases in which one wants to look at a particular charge or $b$-tags, I wrote the pick_charge()and b_jets() functions, which return list-like Objects() classes of leptons of a particular charge and $b$-jets, respectively. E.g.

events[0].b_jets()
events[0]["electron"].pick_charge(1)

This has the advantage that one can still consider list-like objects of e.g. all electrons (positively and negatively charged) and all jets (with and without a $b$-tag). It also matches the LHCO format.

Alternative structure - separate `Objects()` for e.g. lepton and anti-lepton

One alternative is that leptons and anti-leptons, and jets and $b$-jets are treated different objects (i.e. Objects() classes in this case) from the get-go, rather than the same object but with different properties. E.g.

events[0]["b-jet"]
events[0]["anti-electron"]

I'm not a big fan of this, as lots of selections apply for e.g. electrons and anti-electrons. For an $\eta$ cut, you would have to do:

events.cut_objects("electron", eta)
events.cut_objects("anti-electron", eta)  # Got to cut on anti-electrons too

Other possibility

I don't think this is a good idea:

events[0]["electron"]  # All electrons
events[0]["-electron"]  # Only negative charge electrons
events[0]["+electron"]  # Only positive charge electrons

Because inconsistencies might arise when you start removing or counting objects, as objects exist in more than one list. E.g. you might remove electrons with an $\eta$ cut:

events.cut_objects("electron", eta)  # ... but I forgot the "-electron" and "+electron"

and forget about the separate -electron and +electron objects. I don't think this should be possible. The first option gets around this, by defining functions that return e.g. anti-electrons based on the electron objects.

Discussion

Any comments or discussion is welcome.

Error in manual v1 regarding appending events

p14 of v1 of the manual contains a mistake. The offending line of code is:

my_events.append(events[:100]) # Append first 100 events from LHCO file

This ultimately results in an exception: Exception: Cannot overwrite my_events.lhco This should in fact read:

my_events += events[:100]  # Add first 100 events from LHCO file

I will fix this in v2 of the manual.

Bug in $\phi$ function in `Fourvector`

Hi Andrew,

The formula for the phi-function in LHCO_reader here

tan_phi = self[2] / self[1]
phi = np.arctan(tan_phi)

sometimes gives the wrong results. For example,

>>> print(events[0])
+----------+--------+--------+-------+-------+------+------+-------+
|  Object  |  eta   |  phi   |   PT  | jmass | ntrk | btag | hadem |
+----------+--------+--------+-------+-------+------+------+-------+
| electron | -0.501 | 1.769  | 86.42 |  0.0  | -1.0 | 0.0  |  0.0  |
| electron | -0.085 | 2.682  | 51.11 |  0.0  | -1.0 | 0.0  |  0.0  |
| electron | 0.891  | -1.581 | 42.42 |  0.0  | 1.0  | 0.0  |  0.0  |
| electron | 0.681  | -2.481 | 19.46 |  0.0  | 1.0  | 0.0  |  0.0  |
|   muon   | 1.549  | -1.389 | 69.64 |  0.11 | 1.0  | 0.0  |  70.0 |
|   muon   | -0.147 | 0.838  | 11.45 |  0.11 | -1.0 | 0.0  |  11.0 |
|   tau    | 0.983  | -0.283 | 65.58 |  0.0  | 1.0  | 0.0  | 999.9 |
|   tau    | 1.606  | 1.875  | 34.54 |  0.0  | -1.0 | 0.0  | 999.9 |
|   MET    |  0.0   | -0.872 |  8.98 |  0.0  | 0.0  | 0.0  |  0.0  |
+----------+--------+--------+-------+-------+------+------+-------+
>>> e4 = events[0]["electron"][3]
>>> print(e4.vector().phi())
0.66059265359
>>> print(e4.vector()[1],e4.vector()[2])
(-15.366175025038103, -11.940362854616074)

The phi of the 4th electron is -2.481, but phi-function returns 0.66059. The reason is that the values of arctan(tan_phi) are required to be from -0.5*pi to 0.5*pi.

In addition, for the "example.lhco" provided by the LHCO_reader, I found all the phi for the all particles are always in the range 0 to 2*pi.

However, I used the madgraph-pythia-delphes3 to generate LHCO file, in which all the phi are always in the range -pi to pi.

For different convention of the range of phi, phi-function should have different expressions.