diff --git a/.gitattributes b/.gitattributes
new file mode 100644
index 0000000..d9178f8
--- /dev/null
+++ b/.gitattributes
@@ -0,0 +1,9 @@
+*.psd filter=lfs diff=lfs merge=lfs -text
+*.pdf filter=lfs diff=lfs merge=lfs -text
+*.raw filter=lfs diff=lfs merge=lfs -text
+*.root filter=lfs diff=lfs merge=lfs -text
+*.png filter=lfs diff=lfs merge=lfs -text
+*.zip filter=lfs diff=lfs merge=lfs -text
+releases export-ignore
+.gitattributes export-ignore
+.gitignore export-ignore
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..fff1217
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,9 @@
+/GPATH
+/GRTAGS
+/GTAGS
+/__pycache__/
+/old/
+/old1/
+/old2/
+/src/__pycache__/
+src/pygan.egg-info/
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..261eeb9
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,201 @@
+                                 Apache License
+                           Version 2.0, January 2004
+                        http://www.apache.org/licenses/
+
+   TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION
+
+   1. Definitions.
+
+      "License" shall mean the terms and conditions for use, reproduction,
+      and distribution as defined by Sections 1 through 9 of this document.
+
+      "Licensor" shall mean the copyright owner or entity authorized by
+      the copyright owner that is granting the License.
+
+      "Legal Entity" shall mean the union of the acting entity and all
+      other entities that control, are controlled by, or are under common
+      control with that entity. For the purposes of this definition,
+      "control" means (i) the power, direct or indirect, to cause the
+      direction or management of such entity, whether by contract or
+      otherwise, or (ii) ownership of fifty percent (50%) or more of the
+      outstanding shares, or (iii) beneficial ownership of such entity.
+
+      "You" (or "Your") shall mean an individual or Legal Entity
+      exercising permissions granted by this License.
+
+      "Source" form shall mean the preferred form for making modifications,
+      including but not limited to software source code, documentation
+      source, and configuration files.
+
+      "Object" form shall mean any form resulting from mechanical
+      transformation or translation of a Source form, including but
+      not limited to compiled object code, generated documentation,
+      and conversions to other media types.
+
+      "Work" shall mean the work of authorship, whether in Source or
+      Object form, made available under the License, as indicated by a
+      copyright notice that is included in or attached to the work
+      (an example is provided in the Appendix below).
+
+      "Derivative Works" shall mean any work, whether in Source or Object
+      form, that is based on (or derived from) the Work and for which the
+      editorial revisions, annotations, elaborations, or other modifications
+      represent, as a whole, an original work of authorship. For the purposes
+      of this License, Derivative Works shall not include works that remain
+      separable from, or merely link (or bind by name) to the interfaces of,
+      the Work and Derivative Works thereof.
+
+      "Contribution" shall mean any work of authorship, including
+      the original version of the Work and any modifications or additions
+      to that Work or Derivative Works thereof, that is intentionally
+      submitted to Licensor for inclusion in the Work by the copyright owner
+      or by an individual or Legal Entity authorized to submit on behalf of
+      the copyright owner. For the purposes of this definition, "submitted"
+      means any form of electronic, verbal, or written communication sent
+      to the Licensor or its representatives, including but not limited to
+      communication on electronic mailing lists, source code control systems,
+      and issue tracking systems that are managed by, or on behalf of, the
+      Licensor for the purpose of discussing and improving the Work, but
+      excluding communication that is conspicuously marked or otherwise
+      designated in writing by the copyright owner as "Not a Contribution."
+
+      "Contributor" shall mean Licensor and any individual or Legal Entity
+      on behalf of whom a Contribution has been received by Licensor and
+      subsequently incorporated within the Work.
+
+   2. Grant of Copyright License. Subject to the terms and conditions of
+      this License, each Contributor hereby grants to You a perpetual,
+      worldwide, non-exclusive, no-charge, royalty-free, irrevocable
+      copyright license to reproduce, prepare Derivative Works of,
+      publicly display, publicly perform, sublicense, and distribute the
+      Work and such Derivative Works in Source or Object form.
+
+   3. Grant of Patent License. Subject to the terms and conditions of
+      this License, each Contributor hereby grants to You a perpetual,
+      worldwide, non-exclusive, no-charge, royalty-free, irrevocable
+      (except as stated in this section) patent license to make, have made,
+      use, offer to sell, sell, import, and otherwise transfer the Work,
+      where such license applies only to those patent claims licensable
+      by such Contributor that are necessarily infringed by their
+      Contribution(s) alone or by combination of their Contribution(s)
+      with the Work to which such Contribution(s) was submitted. If You
+      institute patent litigation against any entity (including a
+      cross-claim or counterclaim in a lawsuit) alleging that the Work
+      or a Contribution incorporated within the Work constitutes direct
+      or contributory patent infringement, then any patent licenses
+      granted to You under this License for that Work shall terminate
+      as of the date such litigation is filed.
+
+   4. Redistribution. You may reproduce and distribute copies of the
+      Work or Derivative Works thereof in any medium, with or without
+      modifications, and in Source or Object form, provided that You
+      meet the following conditions:
+
+      (a) You must give any other recipients of the Work or
+          Derivative Works a copy of this License; and
+
+      (b) You must cause any modified files to carry prominent notices
+          stating that You changed the files; and
+
+      (c) You must retain, in the Source form of any Derivative Works
+          that You distribute, all copyright, patent, trademark, and
+          attribution notices from the Source form of the Work,
+          excluding those notices that do not pertain to any part of
+          the Derivative Works; and
+
+      (d) If the Work includes a "NOTICE" text file as part of its
+          distribution, then any Derivative Works that You distribute must
+          include a readable copy of the attribution notices contained
+          within such NOTICE file, excluding those notices that do not
+          pertain to any part of the Derivative Works, in at least one
+          of the following places: within a NOTICE text file distributed
+          as part of the Derivative Works; within the Source form or
+          documentation, if provided along with the Derivative Works; or,
+          within a display generated by the Derivative Works, if and
+          wherever such third-party notices normally appear. The contents
+          of the NOTICE file are for informational purposes only and
+          do not modify the License. You may add Your own attribution
+          notices within Derivative Works that You distribute, alongside
+          or as an addendum to the NOTICE text from the Work, provided
+          that such additional attribution notices cannot be construed
+          as modifying the License.
+
+      You may add Your own copyright statement to Your modifications and
+      may provide additional or different license terms and conditions
+      for use, reproduction, or distribution of Your modifications, or
+      for any such Derivative Works as a whole, provided Your use,
+      reproduction, and distribution of the Work otherwise complies with
+      the conditions stated in this License.
+
+   5. Submission of Contributions. Unless You explicitly state otherwise,
+      any Contribution intentionally submitted for inclusion in the Work
+      by You to the Licensor shall be under the terms and conditions of
+      this License, without any additional terms or conditions.
+      Notwithstanding the above, nothing herein shall supersede or modify
+      the terms of any separate license agreement you may have executed
+      with Licensor regarding such Contributions.
+
+   6. Trademarks. This License does not grant permission to use the trade
+      names, trademarks, service marks, or product names of the Licensor,
+      except as required for reasonable and customary use in describing the
+      origin of the Work and reproducing the content of the NOTICE file.
+
+   7. Disclaimer of Warranty. Unless required by applicable law or
+      agreed to in writing, Licensor provides the Work (and each
+      Contributor provides its Contributions) on an "AS IS" BASIS,
+      WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
+      implied, including, without limitation, any warranties or conditions
+      of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A
+      PARTICULAR PURPOSE. You are solely responsible for determining the
+      appropriateness of using or redistributing the Work and assume any
+      risks associated with Your exercise of permissions under this License.
+
+   8. Limitation of Liability. In no event and under no legal theory,
+      whether in tort (including negligence), contract, or otherwise,
+      unless required by applicable law (such as deliberate and grossly
+      negligent acts) or agreed to in writing, shall any Contributor be
+      liable to You for damages, including any direct, indirect, special,
+      incidental, or consequential damages of any character arising as a
+      result of this License or out of the use or inability to use the
+      Work (including but not limited to damages for loss of goodwill,
+      work stoppage, computer failure or malfunction, or any and all
+      other commercial damages or losses), even if such Contributor
+      has been advised of the possibility of such damages.
+
+   9. Accepting Warranty or Additional Liability. While redistributing
+      the Work or Derivative Works thereof, You may choose to offer,
+      and charge a fee for, acceptance of support, warranty, indemnity,
+      or other liability obligations and/or rights consistent with this
+      License. However, in accepting such obligations, You may act only
+      on Your own behalf and on Your sole responsibility, not on behalf
+      of any other Contributor, and only if You agree to indemnify,
+      defend, and hold each Contributor harmless for any liability
+      incurred by, or claims asserted against, such Contributor by reason
+      of your accepting any such warranty or additional liability.
+
+   END OF TERMS AND CONDITIONS
+
+   APPENDIX: How to apply the Apache License to your work.
+
+      To apply the Apache License to your work, attach the following
+      boilerplate notice, with the fields enclosed by brackets "[]"
+      replaced with your own identifying information. (Don't include
+      the brackets!)  The text should be enclosed in the appropriate
+      comment syntax for the file format. We also recommend that a
+      file or class name and description of purpose be included on the
+      same "printed page" as the copyright notice for easier
+      identification within third-party archives.
+
+   Copyright [yyyy] [name of copyright owner]
+
+   Licensed under the Apache License, Version 2.0 (the "License");
+   you may not use this file except in compliance with the License.
+   You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+   Unless required by applicable law or agreed to in writing, software
+   distributed under the License is distributed on an "AS IS" BASIS,
+   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+   See the License for the specific language governing permissions and
+   limitations under the License.
diff --git a/bin/phsp_info b/bin/phsp_info
new file mode 100644
index 0000000..6acef26
--- /dev/null
+++ b/bin/phsp_info
@@ -0,0 +1,51 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import sys
+import numpy as np
+import phsp
+import click
+import os
+
+CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help'])
+@click.command(context_settings=CONTEXT_SETTINGS)
+@click.argument('filename')
+@click.option('-n', default=int(-1), help='Use -1 to read all data')
+def phsp_info(filename, n):
+    '''
+    \b
+    Print information about the given PHSP phase-space file
+
+    \b
+    <FILENAME> : input PHSP root/npy file
+    '''
+
+    data, keys, m = phsp.load(filename, n)
+
+    # print info
+    print('File:       ', filename)
+    b, extension = os.path.splitext(filename)
+    if (extension == '.root'):
+        t = 'root'
+    if (extension == '.npy'):
+        t = 'npy'
+    print('Type:       ', t)
+    print('Nb values:   {0} ({0:.2e})'.format(m))
+    n = len(data)
+    print('Read values: {0} ({0:.2e})'.format(n))
+    print('Nb of keys: ', len(keys))
+    print('Keys:       ', *keys)
+
+    # stats info per key
+    print('{:<10} {:>10} {:>10} {:>10} {:>10}'.format('key', 'min', 'max', 'mean', 'std'))
+    i=0
+    for k in keys:
+        x = data[:,i]
+        print('{:10} {:10.3f} {:10.3f} {:10.3f} {:10.3f}'.format(k, np.amin(x), np.amax(x), np.mean(x), np.std(x)))
+        i = i+1
+
+
+# --------------------------------------------------------------------------
+if __name__ == '__main__':
+    phsp_info()
+
diff --git a/bin/phsp_info~ b/bin/phsp_info~
new file mode 100644
index 0000000..6090d7e
--- /dev/null
+++ b/bin/phsp_info~
@@ -0,0 +1,148 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import sys
+import numpy as np
+import gategan as gg
+#import psf_helpers as h
+import phsp_helpers as h
+import fig_helpers as fh
+import click
+import matplotlib.pyplot as plt
+import os
+from itertools import combinations
+
+CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help'])
+@click.command(context_settings=CONTEXT_SETTINGS)
+@click.argument('filenames', nargs=-1)
+@click.option('--n', '-n',
+              default=1e6,
+              help='Number of samples to generate')
+@click.option('--nb_bins', '-b',
+              default=int(100),
+              help='Number of bins')
+@click.option('--radius', '-r',
+              default=300.0,
+              help='Cylindrical radius in mm')
+@click.option('--mode', '-m', multiple=True,
+              #default='cylinder_1',
+              type=click.Choice(['cylinder_1', 'cylinder_2', 'plane', 'cartesian']),
+              help='Spherical coord parameters mode')
+@click.option('--z_plane', '-z',
+              default=272.1,
+              help='Plane position z in mm')
+def phsp_plot(filenames, nb_bins, n, mode, radius, z_plane):
+    '''
+    \b
+    Compare a TODO
+
+    \b
+    <FILENAME> : input PHSP root/raw file
+
+    '''
+
+    # associate file/mode
+    modes = mode
+    if len(filenames) != len(modes):
+        print('Error, should be as many filename as -m', filenames, modes)
+        exit(0)
+
+    names1 = ['E', 'x', 'dy', 'dz', 'dx']
+    names2 = ['y', 'y', 'y', 'x', 'y']
+
+    # FIXME fig size
+    n = int(n)
+    f, ax = plt.subplots(5, 4, figsize=(15,12))
+
+    # loop on filenames
+    fn=1
+    s=''
+    ranges = {}
+    for filename, mode in zip(filenames,modes):
+        print(filename, mode)
+        value = radius
+        if mode == 'plane' or mode == 'cartesian':
+            value = z_plane
+        print('value', value)
+        f, file_extension = os.path.splitext(filename)
+        names = ['E', 'x', 'y', 'z', 'dx', 'dy', 'dz']
+        if file_extension == '.npy':
+            data_sph, names_sph = h.read_phsp_npy(filename)
+            print(names_sph)
+            m = len(data_sph)
+            data_sph = data_sph[:n]
+            data = h.spherical_to_cartesian(data_sph, value, mode)
+        else:
+            data = h.read_phsp(filename)
+            m = len(data)
+            names_sph = ['E', 'psi', 'z', 'theta', 'phi']
+            if mode == 'plane':
+                names_sph = ['E', 'x', 'y', 'theta', 'phi']
+            if mode == 'cartesian':
+                names_sph = ['E', 'x', 'y', 'dx', 'dy', 'dz']
+            data = data[:n]
+            data_sph = h.cartesian_to_spherical(data, mode)
+        print(names)
+        print(names_sph)
+        print('Number of values:', m, '-> using',n)
+
+        # 1D plot2D
+        i=0
+        for dim in range(len(names)):
+            a = fh.get_sub_fig(ax, i)
+            fh.phsp_plot_histo1D(a, data[:,dim], nb_bins, 'g', str(fn), names[dim])
+            i += 1
+
+        for dim in range(len(names_sph)):
+            na = names_sph[dim]
+            if na in names:
+                print('skip', na)
+            else:
+                a = fh.get_sub_fig(ax, i)
+                fh.phsp_plot_histo1D(a, data_sph[:,dim], nb_bins, 'r', str(fn), names_sph[dim])
+                i += 1
+
+
+        # 2D plot
+        i = i+(len(names1)*(fn-1))
+        for n1, n2 in zip(names1, names2):
+            d1 = data
+            d2 = data
+            i1 = names.index(n1) if n1 in names else None
+            i2 = names.index(n2) if n2 in names else None
+            if i1 == None:
+                i1 = names_sph.index(n1) if n1 in names_sph else None
+                d1 = data_sph
+            if i2 == None:
+                i2 = names_sph.index(n2) if n2 in names_sph else None
+                d2 = data_sph
+            cmap=plt.cm.Reds
+            if i1 == None or i2 == None:
+                print('Error, cannot find ', i1, i2)
+                break
+            a = fh.get_sub_fig(ax, i)
+            nr = n1+n2
+            if nr in ranges:
+                r = ranges[nr]
+            else:
+                r = ''
+            xmin, xmax, ymin, ymax = fh.phsp_plot_histo2D(a, d1, d2, i1, i2, nb_bins, cmap, n1, n2, r)
+            if nr not in ranges:
+                ranges[nr] = [[xmin,xmax],[ymin,ymax]]
+            a.set_title(filename)
+            i += 1
+
+
+        # end loop
+        s += str(fn)+'='+filename+'   '
+        fn = fn+1
+
+    s += '   n='+str(n)
+    plt.suptitle(s)
+    plt.tight_layout()
+    plt.show()
+
+# --------------------------------------------------------------------------
+if __name__ == '__main__':
+    phsp_plot()
+
diff --git a/phsp b/phsp
new file mode 120000
index 0000000..e831038
--- /dev/null
+++ b/phsp
@@ -0,0 +1 @@
+src
\ No newline at end of file
diff --git a/readme.md b/readme.md
new file mode 100644
index 0000000..f8a365e
--- /dev/null
+++ b/readme.md
@@ -0,0 +1,3 @@
+FIXME 
+
+
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000..860d898
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,30 @@
+import setuptools
+
+with open("readme.md", "r") as fh:
+    long_description = fh.read()
+
+setuptools.setup(
+    name="phsp",
+    version="0.0.1",
+    author="David Sarrut",
+    author_email="david.sarrut@creatis.insa-lyon.fr",
+    description="Python tools for GATE phase-space (PHSP) management",
+    long_description=long_description,
+    long_description_content_type="text/markdown",
+    url="https://github.com/dsarrut/phsp",
+    package_dir={'':'src'},
+    packages=setuptools.find_packages('phsp'),
+    classifiers=(
+        "Programming Language :: Python :: 3",
+        "License :: OSI Approved :: Apache Software License",
+        "Operating System :: OS Independent",
+    ),
+    install_requires=[
+        'tqdm',
+        'uproot',
+      ],
+    scripts=[
+        #'bin/phsp_convert',
+        'bin/phsp_info',
+    ]
+)
diff --git a/src/__init__.py b/src/__init__.py
new file mode 100644
index 0000000..b9e5624
--- /dev/null
+++ b/src/__init__.py
@@ -0,0 +1,5 @@
+
+# import files
+from .phsp import *
+from ...bin import *
+
diff --git a/src/__init__.py~ b/src/__init__.py~
new file mode 100644
index 0000000..e697f8f
--- /dev/null
+++ b/src/__init__.py~
@@ -0,0 +1,10 @@
+
+# import files
+from .gategan import *
+from .gategan2 import *
+from .helpers import *
+from .phsp_helpers import *
+from .fig_helpers import *
+from .nn import *
+from ...bin import *
+
diff --git a/src/phsp.egg-info/PKG-INFO b/src/phsp.egg-info/PKG-INFO
new file mode 100644
index 0000000..18652ef
--- /dev/null
+++ b/src/phsp.egg-info/PKG-INFO
@@ -0,0 +1,17 @@
+Metadata-Version: 1.1
+Name: phsp
+Version: 0.0.1
+Summary: Python tools for GATE phase-space (PHSP) management
+Home-page: https://github.com/dsarrut/phsp
+Author: David Sarrut
+Author-email: david.sarrut@creatis.insa-lyon.fr
+License: UNKNOWN
+Description-Content-Type: text/markdown
+Description: FIXME 
+        
+        
+        
+Platform: UNKNOWN
+Classifier: Programming Language :: Python :: 3
+Classifier: License :: OSI Approved :: Apache Software License
+Classifier: Operating System :: OS Independent
diff --git a/src/phsp.egg-info/SOURCES.txt b/src/phsp.egg-info/SOURCES.txt
new file mode 100644
index 0000000..2be614e
--- /dev/null
+++ b/src/phsp.egg-info/SOURCES.txt
@@ -0,0 +1,6 @@
+bin/phsp_info
+src/phsp.egg-info/PKG-INFO
+src/phsp.egg-info/SOURCES.txt
+src/phsp.egg-info/dependency_links.txt
+src/phsp.egg-info/requires.txt
+src/phsp.egg-info/top_level.txt
\ No newline at end of file
diff --git a/src/phsp.egg-info/dependency_links.txt b/src/phsp.egg-info/dependency_links.txt
new file mode 100644
index 0000000..8b13789
--- /dev/null
+++ b/src/phsp.egg-info/dependency_links.txt
@@ -0,0 +1 @@
+
diff --git a/src/phsp.egg-info/requires.txt b/src/phsp.egg-info/requires.txt
new file mode 100644
index 0000000..52a401d
--- /dev/null
+++ b/src/phsp.egg-info/requires.txt
@@ -0,0 +1,2 @@
+tqdm
+uproot
diff --git a/src/phsp.egg-info/top_level.txt b/src/phsp.egg-info/top_level.txt
new file mode 100644
index 0000000..8b13789
--- /dev/null
+++ b/src/phsp.egg-info/top_level.txt
@@ -0,0 +1 @@
+
diff --git a/src/phsp.py b/src/phsp.py
new file mode 100644
index 0000000..0598346
--- /dev/null
+++ b/src/phsp.py
@@ -0,0 +1,99 @@
+#!/bin/env python3
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import sys
+import os
+import uproot
+import time
+
+# -----------------------------------------------------------------------------
+# -----------------------------------------------------------------------------
+# -----------------------------------------------------------------------------
+
+''' ---------------------------------------------------------------------------
+Load a PHSP (Phase-Space) file
+Output is numpy structured array
+'''
+def load(filename, nmax=-1):
+    b, extension = os.path.splitext(filename)
+    nmax = int(nmax)
+
+    if extension == '.root':
+        return load_root(filename, nmax)
+
+    # if extension == '.raw':
+    #     return load_raw(filename)
+
+    if extension == '.npy':
+        return load_npy(filename, nmax)
+
+    print('Error, dont know how to open phsp with extension ',
+          extension,
+          ' (known extensions: .root .npy)')
+    exit(0)
+
+
+''' ---------------------------------------------------------------------------
+Load a PHSP (Phase-Space) file in root format
+Output is numpy structured array
+'''
+def load_root(filename, nmax=-1):
+    nmax = int(nmax)
+    # Check if file exist
+    if (not os.path.isfile(filename)):
+        print("File '"+filename+"' does not exist.")
+        exit()
+
+    # Check if this is a root file
+    try:
+        f = uproot.open(filename)
+    except Exception:
+        print("File '"+filename+"' cannot be opened, not root file ?")
+        exit()
+
+    # Look for a single key named "PhaseSpace"
+    k = f.keys()
+    try:
+        psf = f['PhaseSpace']
+    except Exception:
+        print("This root file does not look like a PhaseSpace, keys are: ",
+              f.keys(), ' while expecting "PhaseSpace"')
+        exit()
+        
+    # Get keys
+    names = [k.decode('UTF-8') for k in psf.keys()]
+    n = psf.numentries
+
+    # Convert to arrays (this take times)
+    if (nmax != -1):
+        a = psf.arrays(entrystop=nmax)
+    else:
+        a = psf.arrays()
+
+    # Concat arrays
+    d = np.column_stack( a[k] for k in psf.keys())
+    
+    return d, names, n
+
+
+''' ---------------------------------------------------------------------------
+Load a PHSP (Phase-Space) file in npy
+Output is numpy structured array
+'''
+def load_npy(filename, nmax=-1):
+    # Check if file exist
+    if (not os.path.isfile(filename)):
+        print("File '"+filename+"' does not exist.")
+        exit()
+
+    x = np.load(filename, mmap_mode='r')
+    n = len(x)
+    print(nmax, n)
+    if nmax > 0:
+        x = x[:nmax]
+    data = x.view(np.float32).reshape(x.shape + (-1,))
+    #data = np.float64(data)
+    return data, x.dtype.names, n
+
+
diff --git a/src/phsp.py~ b/src/phsp.py~
new file mode 100644
index 0000000..031244a
--- /dev/null
+++ b/src/phsp.py~
@@ -0,0 +1,459 @@
+#!/bin/env python3
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import sys
+import os
+import uproot
+from matplotlib import pyplot as plt
+
+# -----------------------------------------------------------------------------
+# -----------------------------------------------------------------------------
+# -----------------------------------------------------------------------------
+
+''' ---------------------------------------------------------------------------
+Read a PHSP (Phase-Space) file
+Output is numpy structured array
+'''
+def read_phsp(filename):
+    b, extension = os.path.splitext(filename)
+
+    if extension == '.root':
+        return read_phsp_root(filename)
+
+    if extension == '.raw':
+        return read_phsp_raw(filename)
+
+    print('Error, dont know how to open phsp with extension ', extension)
+    exit(0)
+
+
+''' ---------------------------------------------------------------------------
+Read a PHSP (Phase-Space) file in root format
+Output is numpy structured array
+'''
+def read_phsp_root(filename):
+   # Check if file exist
+    if (not os.path.isfile(filename)):
+        print("File '"+filename+"' does not exist.")
+        exit()
+
+    # Check if this is a root file
+    try:
+        f = uproot.open(filename)
+    except Exception:
+        print("File '"+filename+"' cannot be opened, not root file ?")
+        exit()
+
+
+    # Look for a single key named "PhaseSpace"
+    k = f.keys()
+    try:
+        psf = f['PhaseSpace']
+    except Exception:
+        print("This root file is not a PhaseSpace, keys are: ", f.keys())
+        exit()
+        psf = f['PhaseSpace']
+
+    # Convert to arrays
+    #print("PhaseSpace keys: ", psf.keys())
+    #psf.show()
+    a = psf.arrays()
+    E = a[b'Ekine']
+    x = a[b'X']
+    y = a[b'Y']
+    z = a[b'Z']
+    dx = a[b'dX']
+    dy = a[b'dY']
+    dz = a[b'dZ']
+
+    # test --> TTree not writable yet
+    # f = uproot.recreate('toto.root')
+    # f['PhaseSpace'] = psf    
+    
+    w = a[b'Weight']
+    uni = np.unique(w)
+    print('Check root weight', uni, len(w), len(uni))
+    if len(uni) > 1:
+        mask = np.where(w == 1.0, 0, 1)
+        E = E[mask==1.0]
+        x = x[mask==1.0]
+        y = y[mask==1.0]
+        z = z[mask==1.0]
+        dx = dx[mask==1.0]
+        dy = dy[mask==1.0]
+        dz = dz[mask==1.0]
+        print('New size is', len(E))
+    
+    data = np.column_stack((E, x, y, z, dx, dy, dz))
+
+    return data
+
+
+''' ---------------------------------------------------------------------------
+Read a PHSP (Phase-Space) file in npy
+Output is numpy structured array
+'''
+def read_phsp_npy(filename):
+   # Check if file exist
+    if (not os.path.isfile(filename)):
+        print("File '"+filename+"' does not exist.")
+        exit()
+
+    x = np.load(filename)
+    data = x.view(np.float32).reshape(x.shape + (-1,))
+    data = np.float64(data)
+    return data, x.dtype.names
+
+
+''' ---------------------------------------------------------------------------
+Read a PHSP (Phase-Space) file in raw
+Output is numpy structured array
+'''
+def read_phsp_raw(filename):
+   # Check if file exist
+    if (not os.path.isfile(filename)):
+        print("File '"+filename+"' does not exist.")
+        exit()
+
+    f = open(filename, "rb")
+    bytes_array = f.read()
+    dt = np.dtype(np.float32)
+    data = np.frombuffer(bytes_array, dtype=dt)
+    rn = len(data)/7
+    data = np.reshape(data, (int(rn),7))
+
+    return data
+
+
+''' ---------------------------------------------------------------------------
+Write data to npy
+'''
+def write_to_npy(data, mode, filename):
+
+    p = []
+    if mode == 'cylinder_1' or mode == 'cylinder_2':
+        r = np.zeros(len(data),
+                     dtype=[('E',     'f4'),
+                            ('psi',   'f4'),
+                            ('z',     'f4'),
+                            ('theta', 'f4'),
+                            ('phi',   'f4'),])
+        r['E'] = data[:,0]
+        r['psi'] = data[:,1]
+        r['z'] = data[:,2]
+        r['theta'] = data[:,3]
+        r['phi'] = data[:,4]
+        np.save(filename, r)
+        return
+
+    if mode == 'plane':
+        r = np.zeros(len(data),
+                     dtype=[('E',     'f4'),
+                            ('x',     'f4'),
+                            ('y',     'f4'),
+                            ('theta', 'f4'),
+                            ('phi',   'f4'),])
+        r['E'] = data[:,0]
+        r['x'] = data[:,1]
+        r['y'] = data[:,2]
+        r['theta'] = data[:,3]
+        r['phi'] = data[:,4]
+        np.save(filename, r)
+        return
+
+    if mode == 'cartesian':
+        r = np.zeros(len(data),
+                     dtype=[('E',   'f4'),
+                            ('x',   'f4'),
+                            ('y',   'f4'),
+                            ('dx',  'f4'),
+                            ('dy',  'f4'),
+                            ('dz',  'f4'),])
+        r['E'] = data[:,0]
+        r['x'] = data[:,1]
+        r['y'] = data[:,2]
+        r['dx'] = data[:,3]
+        r['dy'] = data[:,4]
+        r['dz'] = data[:,5]
+        np.save(filename, r)
+        return
+
+    print('error dont know mode', mode, '(cylinder_1 cylinder_2 plane cartesian)')
+    exit(0)
+
+
+# -----------------------------------------------------------------------------
+# -----------------------------------------------------------------------------
+# -----------------------------------------------------------------------------
+
+''' ---------------------------------------------------------------------------
+Convert to spherical coordinates according to mode
+'''
+def cartesian_to_spherical(data, mode, verbose=False):
+
+    f = { 'cylinder_1': cartesian_to_spherical_cylinder_1,
+          'cylinder_2': cartesian_to_spherical_cylinder_2,
+          'plane': cartesian_to_spherical_plane,
+          'cartesian': cartesian_to_spherical_cartesian
+    }
+
+    if mode in f:
+        return f[mode](data, verbose)
+
+    print('error dont know mode', mode, '(cylinder_1 cylinder_2 plane cartesian)')
+    exit(0)
+
+
+''' ---------------------------------------------------------------------------
+Convert to spherical coordinate according to mode
+'''
+def cartesian_to_spherical_cylinder_1(data, verbose=False):
+
+    E, x, y, z, dx, dy, dz = [data[:,i] for i in np.arange(len(data[0]))]
+
+    radius = np.sqrt(np.power(x,2) + np.power(y, 2))
+    if verbose:
+        print('radius must be constant: ', radius)
+    radius = np.mean(radius)
+
+    psi = np.arctan2(y, x)
+
+    if verbose:
+        r = np.sqrt(dx**2 + dy**2 + dz**2)
+        print('norm dir must be constant: ', r)
+
+    theta = np.arccos(dz)
+    phi = np.arctan2(dy,dx)
+
+    data = np.column_stack((E, psi, z, theta, phi))
+    return data
+
+''' ---------------------------------------------------------------------------
+Convert to spherical coordinate according to mode
+'''
+def cartesian_to_spherical_cylinder_2(data, verbose=False):
+
+    data = cartesian_to_spherical_cylinder_1(data, verbose)
+
+    E, psi, z, theta, phi = [data[:,i] for i in np.arange(len(data))]
+    phi = (phi-psi+np.pi)%(2*np.pi)-np.pi
+
+    data = np.column_stack((E, psi, z, theta, phi))
+    return data
+
+''' ---------------------------------------------------------------------------
+Convert to spherical coordinate according to mode
+'''
+def cartesian_to_spherical_plane(data, verbose=False):
+
+    E, x, y, z, dx, dy, dz = [data[:,i] for i in np.arange(len(data[0]))]
+
+    theta = np.arccos(dz)
+    phi = np.arctan2(dy,dx)
+
+    data = np.column_stack((E, x, y, theta, phi))
+    return data
+
+
+''' ---------------------------------------------------------------------------
+Convert to spherical coordinate according to mode
+'''
+def cartesian_to_spherical_cartesian(data, verbose=False):
+    E, x, y, z, dx, dy, dz = [data[:,i] for i in np.arange(len(data[0]))]
+    data = np.column_stack((E, x, y, dx, dy, dz)) 
+    return data
+
+
+''' ---------------------------------------------------------------------------
+Convert from spherical to cartesian coordinates according to mode
+'''
+def spherical_to_cartesian(data, value, mode, verbose=False):
+
+    f = { 'cylinder_1': spherical_to_cartesian_cylinder_1,
+          'cylinder_2': spherical_to_cartesian_cylinder_2,
+          'plane': spherical_to_cartesian_plane,
+          'cartesian': spherical_to_cartesian_cartesian,}
+
+    if mode in f:
+        return f[mode](data, value, verbose)
+
+    print('error dont know mode', mode, '(cylinder_1 cylinder_2 plane cartesian)')
+    exit(0)
+
+
+''' ---------------------------------------------------------------------------
+Convert from spherical to cartesian coordinates
+'''
+def spherical_to_cartesian_cylinder_1(data, radius, verbose=False):
+
+    E, psi, z, theta, phi = [data[:,i] for i in np.arange(len(data[0]))]
+
+    # position
+    x = radius*np.cos(psi)
+    y = radius*np.sin(psi)
+
+    # convert direction to dx,dy,dz
+    dz = np.cos(theta)
+    f = np.sin(theta)
+    dx = np.cos(phi)
+    dy = np.sin(phi)
+
+    # need to noprm by projected vector lenght
+    # f = np.sin(theta)
+    dy = dy*f
+    dx = dx*f
+
+    data = np.column_stack((E, x, y, z, dx, dy, dz))
+    return data
+
+''' ---------------------------------------------------------------------------
+Convert from spherical to cartesian coordinates
+'''
+def spherical_to_cartesian_cylinder_2(data, radius, verbose=False):
+
+    E, psi, z, theta, phi = [data[:,i] for i in np.arange(len(data[0]))]
+
+    # position
+    x = radius*np.cos(psi)
+    y = radius*np.sin(psi)
+
+    # convert direction to dx,dy,dz
+    phi = (phi+psi+np.pi)%(2*np.pi)-np.pi
+    dz = np.cos(theta)
+    f = np.sin(theta)
+    dx = np.cos(phi)
+    dy = np.sin(phi)
+
+    # need to noprm by projected vector lenght
+    # f = np.sin(theta)
+    dy = dy*f
+    dx = dx*f
+
+    data = np.column_stack((E, x, y, z, dx, dy, dz))
+    return data
+
+''' ---------------------------------------------------------------------------
+Convert from spherical to cartesian coordinates
+'''
+def spherical_to_cartesian_plane(data, z_plane, verbose=False):
+
+    E, x, y, theta, phi = [data[:,i] for i in np.arange(len(data[0]))]
+
+    print('spherical_to_cartesian_plane', z_plane)
+
+    # convert direction to dx,dy,dz
+    dz = np.cos(theta)
+    f = np.sin(theta)
+    dx = np.cos(phi)
+    dy = np.sin(phi)
+
+    # need to noprm by projected vector lenght
+    # f = np.sin(theta)
+    dy = dy*f
+    dx = dx*f
+
+    z = np.ones(len(x))*z_plane
+
+    data = np.column_stack((E, x, y, z, dx, dy, dz))
+    return data
+
+
+''' ---------------------------------------------------------------------------
+Convert from spherical to cartesian coordinates
+'''
+def spherical_to_cartesian_cartesian(data, z_plane, verbose=False):
+    E, x, y, dx, dy, dz = [data[:,i] for i in np.arange(len(data[0]))]
+    z = np.ones(len(x))*z_plane
+    data = np.column_stack((E, x, y, z, dx, dy, dz))
+    return data
+
+
+# -----------------------------------------------------------------------------
+# -----------------------------------------------------------------------------
+# -----------------------------------------------------------------------------
+
+
+''' ---------------------------------------------------------------------------
+Encode all parameters according to the policy
+'''
+def encode_parameters(x, names, policy):
+    ex = np.empty((len(x), 0))
+    param_index = {}
+    for p in policy:
+        a = encode_one_parameter(x, names, p, policy[p])
+        ex = np.column_stack((ex, a))
+        if policy[p] == 'ignore':
+            param_index[p] = -1
+        if policy[p] == 'keep':
+            param_index[p] = len(ex[0])-1
+        if policy[p] == 'sincos':
+            param_index[p] = len(ex[0])-2
+
+    return ex, param_index
+
+
+''' ---------------------------------------------------------------------------
+Decode all parameters according to the policy
+'''
+def decode_parameters(x, names, param_index, policy):
+    ex = np.empty((len(x), 0))
+    for p in names:
+        a = decode_one_parameter(x, param_index[p], p, policy[p])
+        ex = np.column_stack((ex, a))
+
+    return ex
+
+
+''' ---------------------------------------------------------------------------
+Encode one parameter according to the policy
+'''
+def encode_one_parameter(x, names, p, transform):
+    try:
+        index = names.index(p)
+    except:
+        print('error dont know parameter', p)
+        print('in encode_one_parameter', x.shape, names, p, transform)
+        exit(0)
+
+    t = x[:,index]
+
+    if transform == 'ignore':
+        return np.array([])
+
+    if transform == 'keep':
+        return t
+
+    if transform == 'sincos':
+        x1 = np.sin(t)
+        x2 = np.cos(t)
+        t = np.column_stack((x1, x2))
+        return t
+
+    print('error, dont know transform', transform)
+    print('in encode_one_parameter', x.shape, names, p, transform)
+    exit(0)
+
+
+''' ---------------------------------------------------------------------------
+Decode one parameter according to the policy
+'''
+def decode_one_parameter(x, index, p, transform):
+    if transform == 'ignore':
+        return np.array([])
+
+    t = x[:,index]
+
+    if transform == 'keep':
+        return t
+
+    if transform == 'sincos':
+        x1 = t
+        x2 = x[:,index+1]
+        a = np.arctan2(x1,x2)
+        return a
+
+    print('error, dont know transform', transform)
+    print('in decode_one_parameter', x.shape, names, p, transform)
+    exit(0)
+