diff --git a/Control/PyKernel/PyKernel/PyItPatch.h b/Control/PyKernel/PyKernel/PyItPatch.h
new file mode 100755
index 0000000000000000000000000000000000000000..6435c682a7e7d9a5b798aa2ac5c3c11fd7bd2495
--- /dev/null
+++ b/Control/PyKernel/PyKernel/PyItPatch.h
@@ -0,0 +1,74 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef PYANALYSISCORE_PYITPATCH_H
+#define PYANALYSISCORE_PYITPATCH_H
+
+// An utility class to remove '*' from a class type
+template <class T> struct NonPointer
+{
+  typedef T Type;
+};
+template <class T> struct NonPointer<T *>
+{
+  typedef T Type;
+};
+
+// An utility class to remove 'const' from a class type
+template <class T> struct NonConst
+{
+  typedef T Type;
+};
+template <class T> struct NonConst<const T>
+{
+  typedef T Type;
+};
+
+
+/**
+   This class provides some pathces for an iterator. 
+   Some c++ classes have methods to return iterators instead of vectors.
+   In this case people need to use iterators directly. PyKernel addes next()
+   and patches __eq__()/__ne__() instead. Then one can use the iterator like
+   as a python iterator.
+**/
+template <class IT>
+struct PyItPatch
+{
+  PyItPatch() : m_cache(0) {}
+  PyItPatch(IT &it) : m_it(it)
+  { 
+    m_cache = new typename NonConst<typename NonPointer<typename IT::pointer>::Type>::Type();
+  }
+  virtual ~PyItPatch() { if (m_cache!=0) delete m_cache; }
+
+  // next() for python iterator
+  typename IT::reference next()
+  {
+    // this implementation is needed for LCG dict
+    // 'return *m_it++' doesn't compile
+    *m_cache = *m_it;
+    ++m_it;
+    return *m_cache;
+  }
+
+  // __eq__() for python iterator
+  bool eq(IT & rhs)
+  {
+    return rhs == m_it;
+  }
+
+  // __ne__() for python iterator
+  bool ne(IT & rhs)
+  {
+    return !eq(rhs);
+  }
+
+private:
+  IT m_it;
+  typename NonConst<typename NonPointer<typename IT::pointer>::Type>::Type *m_cache;
+};
+
+#endif
+
diff --git a/Control/PyKernel/PyKernel/PyKernelDict.h b/Control/PyKernel/PyKernel/PyKernelDict.h
new file mode 100755
index 0000000000000000000000000000000000000000..efc2b19c0d2de49c8b0cc1d386ffa6801ccb921b
--- /dev/null
+++ b/Control/PyKernel/PyKernel/PyKernelDict.h
@@ -0,0 +1,10 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef PYKERNEL_PYKERNELDICT_H
+#define PYKERNEL_PYKERNELDICT_H
+
+#include "PyKernel/PyReverseProxy.h"
+
+#endif
diff --git a/Control/PyKernel/PyKernel/PyReverseProxy.h b/Control/PyKernel/PyKernel/PyReverseProxy.h
new file mode 100755
index 0000000000000000000000000000000000000000..1795bb10482e6ce9e1ea11c5bba185d0fc694930
--- /dev/null
+++ b/Control/PyKernel/PyKernel/PyReverseProxy.h
@@ -0,0 +1,80 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef PYKERNEL_PYREVERSEPROXY_H
+#define PYKERNEL_PYREVERSEPROXY_H
+
+/**
+   A utility class to convert a C++ object to a PyObject
+
+   @author Tadashi Maeno
+*/
+
+#include "Python.h"
+#include <map>
+#include <string>
+
+
+struct PyReverseProxy
+{
+  /// destructor
+  ~PyReverseProxy()
+  {
+    m_proxyMap.erase(m_key);
+  }
+
+  /// setter
+  void setFunc (PyObject * func)
+  {
+    Py_DECREF(m_func); 
+    m_func = func;
+    Py_INCREF(m_func); 
+  }
+
+  // getter
+  void getObj (void *& obj)
+  {
+    obj = m_obj;
+  }
+
+  PyObject * toPyObj (void * obj)
+  {
+    m_obj = obj;
+    return PyObject_CallObject(m_func,NULL);
+  }
+
+  /// factory method
+  static PyReverseProxy * getProxy(const std::string & key)
+  {
+    if (m_proxyMap.count(key) == 0)
+      m_proxyMap[key] = new PyReverseProxy(key);
+    return m_proxyMap[key];
+  }
+
+private:
+  /// default constructor : never used
+  PyReverseProxy() 
+  {}
+
+  /// constructor
+  PyReverseProxy(const std::string &key)
+    : m_key(key), m_obj(0), m_func(0)
+  {}
+
+  /// key
+  std::string m_key;
+
+  /// C++ obj
+  void * m_obj;
+
+  /// python code fragment to convert C++ obj to PyObj
+  PyObject * m_func;
+
+  /// proxy map
+  static std::map<std::string, PyReverseProxy *> m_proxyMap;
+};
+
+#endif
+
+  
diff --git a/Control/PyKernel/PyKernel/selection.xml b/Control/PyKernel/PyKernel/selection.xml
new file mode 100755
index 0000000000000000000000000000000000000000..98eadef58eb67206d3afd109e086482a4b251ddd
--- /dev/null
+++ b/Control/PyKernel/PyKernel/selection.xml
@@ -0,0 +1,5 @@
+<lcgdict>
+
+  <class name="PyReverseProxy" />
+  
+</lcgdict>
\ No newline at end of file
diff --git a/Control/PyKernel/cmt/requirements b/Control/PyKernel/cmt/requirements
new file mode 100755
index 0000000000000000000000000000000000000000..34185b6364f00150fc825e5f7b80942b3925e4f2
--- /dev/null
+++ b/Control/PyKernel/cmt/requirements
@@ -0,0 +1,19 @@
+package PyKernel
+
+author Tadashi Maeno <Tadashi.Maeno@cern.ch>
+
+use AtlasPolicy		AtlasPolicy-*
+use AtlasPython		AtlasPython-*	External
+
+apply_pattern declare_python_modules files="*.py"
+apply_pattern declare_joboptions     files="*.py"
+
+library PyKernel *.cxx
+apply_pattern installed_library
+
+private
+
+use AtlasReflex  AtlasReflex*  External
+
+apply_pattern lcgdict dict=PyKernel selectionfile=selection.xml \
+        headerfiles="../PyKernel/PyKernelDict.h"
diff --git a/Control/PyKernel/doc/mainpage.h b/Control/PyKernel/doc/mainpage.h
new file mode 100755
index 0000000000000000000000000000000000000000..ac0cf1c3fc7af95ef582584c8c3897192c26fff7
--- /dev/null
+++ b/Control/PyKernel/doc/mainpage.h
@@ -0,0 +1,43 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/**
+
+@mainpage
+
+@section PyKernelIntro Introduction 
+
+The PyKernel package contains some python modules 
+and python-bindings for interactive anaysis using Athena.
+Note that this package is for python but doxygen doesn't
+treat python very well. See <a href="http://tmaeno.home.cern.ch/tmaeno/PyKernel/index.html">PyKernel doc</a> instead.
+
+
+@section PyKernelPythonMod Python modules
+
+  - PyKernel.py: core module to help StoreGate access, histogramming and visualization
+  - PyKHist.py: wrapper module for histograms
+
+@section PyKernelPythonBind Python bindings
+
+  - PyItPatch: contains patches to treat C++ iterators as python's iterator
+  - PyReverseProxy: converts a C++ object to a PyObject
+
+@section PyKernelExtras Extra Pages
+
+  - @ref PyKernelUsed
+  - @ref PyKernelReq
+*/
+
+/**
+@page PyKernelUsed Used Packages
+@htmlinclude used_packages.html
+*/
+
+
+/**
+@page PyKernelReq Requirements
+@include requirements
+
+*/
diff --git a/Control/PyKernel/python/PyKHist.py b/Control/PyKernel/python/PyKHist.py
new file mode 100755
index 0000000000000000000000000000000000000000..648e36f6c9be517ff1f32f6f4327cecd0cb1f462
--- /dev/null
+++ b/Control/PyKernel/python/PyKHist.py
@@ -0,0 +1,198 @@
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+import re
+import math
+import ROOT
+
+# a function object for wrapping ROOT methods
+#     update ROOT histogram and then call the ROOT method
+#   :param pykObj    : PyKHist object  
+#   :param rootClass : base class of the object
+#   :param methodName: method name (e.g., Draw)
+#   :param var       : arguments of the method
+class _PyKMethod:
+    def __init__(self,pykObj,rootClass,methodName):
+        self.pykObj = pykObj
+        # get pointer to the method
+        self.methodObj = getattr(rootClass,methodName)
+    # method emulation
+    def __call__(self,*var):
+        # update ROOT histogram
+        self.pykObj._updateROOT()
+        # call ROOT method
+        return apply(self.methodObj,(self.pykObj,)+var)
+
+
+# special methods not to be patched
+_specialMethod = ['Fill','Add']
+
+
+# wrapper class for 1D histogram
+class PyKHist1D (ROOT.TH1D):
+
+    # constructor
+    #   :param aidaHist: AIDA hist obj
+    def __init__(self,aidaHist):
+        self.aidaHist = aidaHist
+        # go to the global directory. Otherwise this object is destroyed
+        #   when a pool is closed
+        ROOT.gROOT.cd()
+        # init ROOT hist
+        ROOT.TH1D.__init__(self,aidaHist.title(),aidaHist.title(),\
+                           aidaHist.axis().bins(),aidaHist.axis().lowerEdge(),\
+                           aidaHist.axis().upperEdge())
+        # update ROOT histogram
+        self._updateROOT()
+        # patch ROOT methods
+        rootClass = ROOT.TH1D
+        attrs = dir(rootClass)
+        for attr in attrs:
+            # choose ROOT method only
+            if re.match('^[A-Z]',attr) == None:
+                continue
+            attrObj = getattr(rootClass,attr)
+            # method?
+            if not callable(attrObj):
+                continue
+            # special method?
+            if _specialMethod.count(attr) != 0:
+                continue
+            # instantiate function object
+            mObj = _PyKMethod(self,rootClass,attr)
+            # patch
+            setattr(self,attr,mObj)
+
+    # update ROOT part
+    def _updateROOT(self):
+        for i in range(0,self.aidaHist.axis().bins()):
+            ROOT.TH1D.SetBinContent(self,i+1,self.aidaHist.binHeight(i))
+        ROOT.TH1D.SetEntries(self,self.aidaHist.allEntries())
+        
+    # overload of Fill    
+    def Fill(self,*var):
+        apply(self.aidaHist.fill,var)
+        apply(ROOT.TH1D.Fill,(self,)+var)
+
+    # overload of Add
+    def Add(self,hist):
+        self.aidaHist.add(hist.aidaHist)
+        # ROOT.TH1D.Add(self,hist,1) causes a crash       
+        ROOT.TH1D.Add(self,hist,hist,1,0)
+
+
+# wrapper class for 2D histogram
+class PyKHist2D (ROOT.TH2D):
+
+    # constructor
+    #   :param aidaHist: AIDA hist obj
+    def __init__(self,aidaHist):
+        self.aidaHist = aidaHist
+        # go to the global directory. Otherwise this object is destroyed
+        #   when a pool is closed
+        ROOT.gROOT.cd()
+        # init ROOT hist
+        ROOT.TH2D.__init__(self,aidaHist.title(),aidaHist.title(),\
+                           aidaHist.xAxis().bins(),aidaHist.xAxis().lowerEdge(),\
+                           aidaHist.xAxis().upperEdge(),aidaHist.yAxis().bins(),\
+                           aidaHist.yAxis().lowerEdge(),aidaHist.yAxis().upperEdge())
+        # update ROOT histogram
+        self._updateROOT()
+        # patch ROOT methods
+        rootClass = ROOT.TH2D
+        attrs = dir(rootClass)
+        for attr in attrs:
+            # choose ROOT method only
+            if re.match('^[A-Z]',attr) == None:
+                continue
+            attrObj = getattr(rootClass,attr)
+            # method?
+            if not callable(attrObj):
+                continue
+            # special method?
+            if _specialMethod.count(attr) != 0:
+                continue
+            # instantiate function object
+            mObj = _PyKMethod(self,rootClass,attr)
+            # patch
+            setattr(self,attr,mObj)
+                                                                  
+    # update ROOT part
+    def _updateROOT(self):
+        for i in range(0,self.aidaHist.xAxis().bins()):
+            for j in range(0,self.aidaHist.yAxis().bins()):
+                ROOT.TH2D.SetBinContent(self,i+1,j+1,self.aidaHist.binHeight(i,j))
+        ROOT.TH2D.SetEntries(self,self.aidaHist.allEntries())
+
+    # overload of Fill    
+    def Fill(self,*var):
+        apply(self.aidaHist.fill,var)
+        apply(ROOT.TH2D.Fill,(self,)+var)
+
+    # overload of Add
+    def Add(self,hist):
+        self.aidaHist.add(hist.aidaHist)
+        # ROOT.TH2D.Add(self,hist,1) causes a crash       
+        ROOT.TH2D.Add(self,hist,hist,1,0)
+
+
+# wrapper class for profile histogram
+class PyKHistProf (ROOT.TProfile):
+
+    # constructor
+    #   :param aidaHist: AIDA hist obj
+    def __init__(self,aidaHist):
+        self.aidaHist = aidaHist
+        # go to the global directory. Otherwise this object is destroyed
+        #   when a pool is closed
+        ROOT.gROOT.cd()
+        # init ROOT hist
+        ROOT.TProfile.__init__(self,aidaHist.title(),aidaHist.title(),\
+                               aidaHist.axis().bins(),aidaHist.axis().lowerEdge(),\
+                               aidaHist.axis().upperEdge())
+        # update ROOT histogram
+        self._updateROOT()
+        # patch ROOT methods
+        rootClass = ROOT.TProfile
+        attrs = dir(rootClass)
+        for attr in attrs:
+            # choose ROOT method only
+            if re.match('^[A-Z]',attr) == None:
+                continue
+            attrObj = getattr(rootClass,attr)
+            # method?
+            if not callable(attrObj):
+                continue
+            # public method?
+            if attr[0] == '_':
+                continue
+            # special method?
+            if _specialMethod.count(attr) != 0:
+                continue
+            # instantiate function object
+            mObj = _PyKMethod(self,rootClass,attr)
+            # patch
+            setattr(self,attr,mObj)
+
+    # update ROOT part
+    def _updateROOT(self):
+        for i in range(0,self.aidaHist.axis().bins()):
+            sumwyBin = self.aidaHist.binHeight(i)*self.aidaHist.binEntries(i)
+            sumwy2Bin = (self.aidaHist.binRms(i)*self.aidaHist.binRms(i)+\
+                         self.aidaHist.binHeight(i)*self.aidaHist.binHeight(i))*self.aidaHist.binEntries(i)
+            ROOT.TProfile.SetBinContent(self,i+1,sumwyBin)
+            ROOT.TProfile.SetBinError(self,i+1,math.sqrt(sumwy2Bin))
+            ROOT.TProfile.SetBinEntries(self,i+1,self.aidaHist.binEntries(i))
+        ROOT.TProfile.SetEntries(self,self.aidaHist.allEntries())
+                                                      
+    # overload of Fill    
+    def Fill(self,*var):
+        apply(self.aidaHist.fill,var)
+        apply(ROOT.TProfile.Fill,(self,)+var)
+                                            
+    # overload of Add
+    def Add(self,hist):
+        self.aidaHist.add(hist.aidaHist)
+        # ROOT.TProfile.Add(self,hist,1) causes a crash       
+        ROOT.TProfile.Add(self,hist,hist,1,0)
+
+
diff --git a/Control/PyKernel/python/PyKernel.py b/Control/PyKernel/python/PyKernel.py
new file mode 100755
index 0000000000000000000000000000000000000000..2af5df8e0215d9523087269efcebd69dd1b42ba2
--- /dev/null
+++ b/Control/PyKernel/python/PyKernel.py
@@ -0,0 +1,915 @@
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+"""core module for an interactive analysis
+
+Examples are in `PyAnalysisExamples/PlotTest.py`_ and `PyAnalysisExamples/HistTest.py`_
+
+.. _PyAnalysisExamples/HistTest.py: http://atlas-sw.cern.ch/cgi-bin/viewcvs-atlas.cgi/offline/PhysicsAnalysis/PyAnalysis/PyAnalysisExamples/share/HistTest.py?rev=HEAD&content-type=text/vnd.viewcvs-markup
+
+.. _PyAnalysisExamples/PlotTest.py: http://atlas-sw.cern.ch/cgi-bin/viewcvs-atlas.cgi/offline/PhysicsAnalysis/PyAnalysis/PyAnalysisExamples/share/PlotTest.py?rev=HEAD&content-type=text/vnd.viewcvs-markup
+    
+:author: Tadashi Maeno
+:contact: tmaeno@bnl.gov
+
+"""
+__docformat__ = "restructuredtext en"
+
+# remove these lines when epydoc
+#"""
+import re
+import types
+import PyCintex
+import PyKHist
+from math import *
+from AthenaCommon.SystemOfUnits import *
+
+# global name space
+GNS = PyCintex.Namespace('')
+
+#"""
+
+# Event Loop Types
+_PreProcess    = 0
+_NormalProcess = 1
+_HybridProcess = 2
+
+# dummy class for TAG and AANT
+class _DummyClass: pass
+setattr (GNS, 'AttributeList', _DummyClass)
+setattr (GNS, 'AANT',          _DummyClass)
+del _DummyClass
+
+storeGate = None
+detStore = None
+
+
+# initialize core
+def init (v_theApp, v_rootStream=None):
+    """Initialize core
+
+    This method is called in `PyKernel/InitPyKernel.py`_.
+    
+    :param v_theApp: reference to the application manager. theApp
+
+    **examples**::
+
+      athena> PyKernel.init(theApp)
+
+    .. _PyKernel/InitPyKernel.py: http://atlas-sw.cern.ch/cgi-bin/viewcvs-atlas.cgi/offline/Control/PyKernel/share/InitPyKernel.py?rev=HEAD&content-type=text/vnd.viewcvs-markup  
+      
+    """
+    # application manager
+    global theApp
+    theApp = v_theApp
+    # root stream
+    global rootStream
+    rootStream = v_rootStream
+    # event loop type
+    global eventLoopType
+    eventLoopType = _HybridProcess
+    # event counter
+    global curEvent
+    curEvent = 0
+    if hasattr (theApp,  'curEvent'):
+        # patch some methods
+        method = "nextEvent"
+        mobj = _SetEventCounter(method)
+        setattr(theApp.__class__,method,mobj)
+        method = "run"
+        mobj = _SetEventCounter(method)
+        setattr(theApp.__class__,method,mobj)
+
+
+# a function object to set event counter in ROOT stream
+class _SetEventCounter:
+    """Function object to set event counter in ROOT stream.
+    This class is used to patch, e.g., theApp.nextEvent(), run() ...
+
+    :param methodObj: reference to method to be patched
+
+    """
+    def __init__(self,methodName):
+        self.methodObj   = getattr(theApp,methodName)
+    # method emulation    
+    def __call__ (self,*var):
+        global curEvent
+        # set event counter. For pre-process internal curEvent is used
+        if eventLoopType != _PreProcess:
+            curEvent = theApp.curEvent()
+        # get ROOT entry    
+        if eventLoopType != _NormalProcess and rootStream:
+            rootStream.GetEntry(curEvent)
+        # return if pre-process
+        if eventLoopType == _PreProcess:
+            # increment event counter
+            curEvent = curEvent+1
+            return GNS.StatusCode(1)
+        # call original method
+        return apply(self.methodObj,var)
+    
+
+# retrieve object from StoreGate
+def retrieve (aClass, aKey=None):
+    """Retrieve object from StoreGate
+    
+    :param aClass: type of the class
+    :param aKey: key of the object
+
+    **examples**::
+
+      athena> obj = PyKernel.retrieve(g.MyClass,'mykey')
+      athena> obj = PyKernel.retrieve(g.MyClass)  # when only one MyClass obj is in SG
+          where the prefix 'g' is the global namespace provided by PyCintex
+          g = PyCintex.Namespace('')
+
+    """
+    #import workaround
+    if aClass == GNS.AttributeList:
+        return rootStream        
+    if aClass == GNS.AANT:
+        return rootStream
+    global storeGate
+    if storeGate == None:
+        import StoreGateBindings.Bindings
+        storeGate = StoreGateBindings.Bindings.StoreGate.pointer()
+    if aKey:
+        ret = storeGate.retrieve(aClass,aKey)
+    else:
+        ret = storeGate.retrieve(aClass)
+    return ret
+
+
+# retrieve object from DetectorStore
+def retrieveDet (aClass, aKey=None):
+    """Retrieve object from DetectorStore
+    
+    :param aClass: type of the class
+    :param aKey: key of the object
+
+    **examples**::
+    
+      athena> obj = PyKernel.retrieveDet(g.MyClass,'mykey')
+      athena> obj = PyKernel.retrieveDet(g.MyClass) # when only one MyClass obj is in SG
+          where the prefix 'g' is the global namespace provided by PyCintex
+          g = PyCintex.Namespace('')
+          
+    """
+    #import workaround    
+    global detStore
+    if detStore == None:
+        import StoreGateBindings.Bindings
+        detStore = StoreGateBindings.Bindings.StoreGate.pointer("DetectorStore")
+    if aKey:
+        ret = detStore.retrieve(aClass,aKey)
+    else:
+        ret = detStore.retrieve(aClass)        
+    return ret
+
+
+# fill a histogram
+def fill (hist, classAndKey, value, criteria="True", nEvent=100):
+    '''
+    Fill 1D-histogram
+    
+    :param hist: reference to AIDA or ROOT histogram
+    :param classAndKey: combination of class name and key separeted with "#". "Class#Key"
+    :param value: physics parameter in string
+    :param criteria: selection criteria
+    :param nEvent: number of event to be processed
+
+    **examples**::
+
+      athena> fill(h,"ElectronContainer#ElectronCollection","$x.pt()")
+          fill hist with pt of electrons
+         "$x" denotes an element of "Class#Key", if "Class#Key" is a collection
+       
+      athena> fill(h,"MissingET#MET_Calib","$x.sumet()")
+          fill hist with et of MissingET.
+          "$x" denotes "Class#Key" itself, if "Class#Key" is not a vector-like class
+          
+      athena> fill(h,"ElectronContainer#ElectronCollection","$x.pt()","$x.pz()>0")
+          apply a selection criteria
+          
+    For more detail of parameters, see `PyAnalysisExamples/PlotTest.py`_
+
+    .. _PyAnalysisExamples/PlotTest.py: http://atlas-sw.cern.ch/cgi-bin/viewcvs-atlas.cgi/offline/PhysicsAnalysis/PyAnalysis/PyAnalysisExamples/share/PlotTest.py?rev=HEAD&content-type=text/vnd.viewcvs-markup
+    
+    '''
+
+    # number of buffered events
+    bufEvent = nEvent
+    if nEvent > 100:
+        bufEvent = 100
+
+    # convert class&key to a store gate access
+    commandSG = "obj = None"
+    if classAndKey != "":
+        commandSG = "obj = "+_parseString(classAndKey)
+
+    # build commands
+    if callable(value):
+        # if function pointer
+        commandV = "vX = value()"
+    else:
+        # if string, parse value/criteria     
+        commandV = "vX = "+_parseString(value)
+    if callable(criteria):
+        # if function pointer        
+        commandC = "vC = criteria()"
+    else:
+        # if string, parse value/criteria             
+        commandC = "vC = "+_parseString(criteria)
+
+    # initialize application mgr
+    theApp.initialize()
+
+    # buffer to determine x-range of histgram
+    buf = []
+
+    # loop over nEvent
+    for iE in range(nEvent):
+        # get object from SG
+        theApp.nextEvent()
+        try:
+            exec commandSG
+
+            # check if the obj is a vector-like class
+            if hasattr(obj,'size') and hasattr(obj,'__getitem__'):
+                lSize = obj.size()
+                isCollection = True
+            else:
+                lSize = 1
+                isCollection = False
+            # if NULL
+            if obj == 0:
+                lSize = 0                
+        except:
+            lSize = 0
+
+        # loop over all elements
+        for iC in range(lSize):
+            # parameter name "x" must be consistent with the parsed strings
+            if isCollection:
+                x = obj[iC]
+            else:
+                x = obj
+
+            # exec value/criteria commands
+            try:
+                exec commandV
+                exec commandC
+
+                # evaluate vC/vX. "vC" and "vX" must be consistent with commands
+                if vC:
+                    if iE < bufEvent:
+                        buf.append(vX)
+                    else:
+                        h.Fill(vX)
+            except:
+                pass
+
+            # if not a collection escape from loop
+            if not isCollection:
+                break
+            
+        # create Histogram
+        if (iE+1) == bufEvent:
+            if len(buf)==0:
+                minX=0
+            else:    
+                minX = min(buf)
+                
+            if minX<0:
+                minX *= 1.1
+            elif minX>0:
+                minX *= 0.9
+            else:
+                minX = -1
+
+            if len(buf)==0:
+                maxX=0
+            else:    
+                maxX = max(buf)
+                
+            if maxX<0:
+                maxX *= 0.9
+            elif maxX>0:
+                maxX *= 1.1
+            else:
+                maxX = 1
+
+            # create histogram if hist is None
+            if hist == None:
+                lpath = '/stat/tmp/PyKernelHist'
+                unregister(lpath)
+                # determine title of histo
+                if callable(value):
+                    title = value.__name__
+                else:
+                    title = value
+                h = book(lpath, title, 100, minX, maxX)
+            else:
+                h = hist
+
+    # buffered elements
+    for vB in buf:
+        h.Fill(vB)
+
+    return h
+
+
+# plot a histogram
+def plot (classAndKey, value="$x", criteria="True", nEvent=100):
+    """Plot 1D-histogram
+    
+    :param classAndKey: combination of class name and key separeted with '#'. 'Class#Key'
+    :param value: physics parameter in string
+    :param criteria: selection criteria
+    :param nEvent: number of event to be processed
+
+    **examples**::
+
+      athena> plot('ElectronContainer#ElectronCollection','$x.pt()')
+          plot pt of electrons
+          
+    For detail, see `PyAnalysisExamples/PlotTest.py`_
+
+    .. _PyAnalysisExamples/PlotTest.py: http://atlas-sw.cern.ch/cgi-bin/viewcvs-atlas.cgi/offline/PhysicsAnalysis/PyAnalysis/PyAnalysisExamples/share/PlotTest.py?rev=HEAD&content-type=text/vnd.viewcvs-markup
+    
+    """
+  
+    # fill a histogram
+    h = fill (None, classAndKey, value, criteria, nEvent)
+    # draw 
+    h.Draw()
+    # this return is needed to draw up the canvas.
+    # note that PyKernel._hSave = h doesn't work although it keeps the hist in memory  
+    return h
+
+
+# fill 2D histogram
+def fill2 (hist, classAndKey, valueX, valueY, criteria="True", nEvent=100):
+    """Fill 2D-histogram
+    
+    :param hist: reference to AIDA or ROOT histogram
+    :param classAndKey: combination of class name and key separeted with '#', 'Class#Key'
+    :param valueX: physics parameter for X in string
+    :param valueY: physics parameter for Y in string
+    :param criteria: selection criteria
+    :param nEvent: number of event to be processed
+
+    For detail, see `fill`
+    
+    """
+
+    # number of buffered events
+    bufEvent = nEvent
+    if nEvent > 100:
+        bufEvent = 100
+
+    # convert class&key to a store gate access
+    commandSG = "obj = None"
+    if classAndKey != "":
+        commandSG = "obj = "+_parseString(classAndKey)
+    
+    # build commands
+    if callable(valueX):
+        # if function pointer
+        commandX = "vX = valueX()"
+    else:
+        # if string, parse value/criteria     
+        commandX = "vX = "+_parseString(valueX)
+    if callable(valueY):
+        # if function pointer
+        commandY = "vY = valueY()"
+    else:
+        # if string, parse value/criteria     
+        commandY = "vY = "+_parseString(valueY)
+    if callable(criteria):
+        # if function pointer        
+        commandC = "vC = criteria()"
+    else:
+        # if string, parse value/criteria             
+        commandC = "vC = "+_parseString(criteria)
+
+    # initialize application mgr
+    theApp.initialize()
+
+    # buffer to determine xy-range of histgram
+    bufX = []
+    bufY = []    
+    
+    # loop over nEvent
+    for iE in range(nEvent):
+        # get object from SG
+        theApp.nextEvent()        
+        try:
+            exec commandSG
+
+            # check if the obj is a vector-like class
+            if hasattr(obj,'size') and hasattr(obj,'__getitem__'):
+                lSize = obj.size()
+                isCollection = True
+            else:
+                lSize = 1
+                isCollection = False
+            # if NULL
+            if obj == 0:
+                lSize = 0                
+        except:
+            lSize = 0            
+        
+        # loop over all elements
+        for iC in range(lSize):
+            # parameter name "x" must be consistent with the parsed strings
+            if isCollection:
+                x = obj[iC]
+            else:
+                x = obj
+
+            # exec value/criteria commands
+            try:
+                exec commandX
+                exec commandY                
+                exec commandC
+
+                # evaluate vC/vX. "vC" and "vX" must be consistent with commands
+                if vC:
+                    if iE < bufEvent:
+                        bufX.append(vX)
+                        bufY.append(vY)                        
+                    else:
+                        h.Fill(vX,vY)
+            except:
+                pass
+                            
+            # if not a collection escape from loop
+            if not isCollection:
+                break
+
+        # create Histogram
+        if (iE+1) == bufEvent:
+            if len(bufX)==0:
+                minX=0
+            else:    
+                minX = min(bufX)
+                
+            if minX<0:
+                minX *= 1.1
+            elif minX>0:
+                minX *= 0.9
+            else:
+                minX = -1
+                
+            if len(bufX)==0:
+                maxX=0
+            else:    
+                maxX = max(bufX)
+                
+            if maxX<0:
+                maxX *= 0.9
+            elif maxX>0:
+                maxX *= 1.1
+            else:
+                maxX = 1
+
+            if len(bufY)==0:
+                minY=0
+            else:    
+                minY = min(bufY)
+                
+            if minY<0:
+                minY *= 1.1
+            elif minY>0:
+                minY *= 0.9
+            else:
+                minY = -1
+
+            if len(bufY)==0:
+                maxY=0
+            else:    
+                maxY = max(bufY)
+                
+            if maxY<0:
+                maxY *= 0.9
+            elif maxY>0:
+                maxY *= 1.1
+            else:
+                maxY = 1
+
+            # create histogram if hist is None
+            if hist == None:
+                lpath = '/stat/tmp/PyKernelHist'
+                unregister(lpath)
+                # determine title of histo
+                if callable(valueX):
+                    titleX = valueX.__name__
+                else:
+                    titleX = valueX                    
+                if callable(valueY):
+                    titleY = valueY.__name__
+                else:
+                    titleY = valueY
+                h = book(lpath, titleY+" vs "+titleX, 100, minX, maxX, 100, minY, maxY)
+            else:
+                h = hist
+
+    # buffered elements
+    for iB in range(len(bufX)):
+        vX = bufX[iB]
+        vY = bufY[iB]
+        h.Fill(vX,vY)
+        
+    return h    
+
+
+# plot 2D histogram
+def plot2 (classAndKey, valueX="$x", valueY="$x", criteria="True", nEvent=100):
+    """Plot 2D-histogram
+    
+    :param classAndKey: combination of class name and key separeted with '#', 'Class#Key'
+    :param valueX: physics parameter for X in string
+    :param valueY: physics parameter for Y in string
+    :param criteria: selection criteria
+    :param nEvent: number of event to be processed
+
+    For detail, see `plot`
+    
+    """
+    h = fill2 (None, classAndKey, valueX, valueY, criteria, nEvent)
+    # draw 
+    h.Draw('BOX')
+    # this return is needed to draw up the canvas.
+    # note that PyKernel._hSave = h doesn't work although it keeps the hist in memory  
+    return h
+
+    
+# fill profile histogram
+def fillProf (hist, classAndKey, valueX, valueY, criteria="True", nEvent=100):
+    """Fill profile-histogram
+    
+    :param hist: reference to AIDA or ROOT histogram
+    :param classAndKey: combination of class name and key separeted with '#', 'Class#Key'
+    :param valueX: physics parameter for X in string
+    :param valueY: physics parameter for Y in string
+    :param criteria: selection criteria
+    :param nEvent: number of event to be processed
+
+    For detail, see `fill`
+    
+    """
+
+    # number of buffered events
+    bufEvent = nEvent
+    if nEvent > 100:
+        bufEvent = 100
+
+    # convert class&key to a store gate access
+    commandSG = "obj = None"
+    if classAndKey != "":
+        commandSG = "obj = "+_parseString(classAndKey)
+    
+    # build commands
+    if callable(valueX):
+        # if function pointer
+        commandX = "vX = valueX()"
+    else:
+        # if string, parse value/criteria     
+        commandX = "vX = "+_parseString(valueX)
+    if callable(valueY):
+        # if function pointer
+        commandY = "vY = valueY()"
+    else:
+        # if string, parse value/criteria     
+        commandY = "vY = "+_parseString(valueY)
+    if callable(criteria):
+        # if function pointer        
+        commandC = "vC = criteria()"
+    else:
+        # if string, parse value/criteria             
+        commandC = "vC = "+_parseString(criteria)
+
+    # initialize application mgr
+    theApp.initialize()
+
+    # buffer to determine xy-range of histgram
+    bufX = []
+    bufY = []    
+    
+    # loop over nEvent
+    for iE in range(nEvent):
+        # get object from SG
+        theApp.nextEvent()        
+        try:
+            exec commandSG
+
+            # check if the obj is a vector-like class
+            if hasattr(obj,'size') and hasattr(obj,'__getitem__'):
+                lSize = obj.size()
+                isCollection = True
+            else:
+                lSize = 1
+                isCollection = False
+            # if NULL
+            if obj == 0:
+                lSize = 0                
+        except:
+            lSize = 0            
+        
+        # loop over all elements
+        for iC in range(lSize):
+            # parameter name "x" must be consistent with the parsed strings
+            if isCollection:
+                x = obj[iC]
+            else:
+                x = obj
+
+            # exec value/criteria commands
+            try:
+                exec commandX
+                exec commandY                
+                exec commandC
+
+                # evaluate vC/vX. "vC" and "vX" must be consistent with commands
+                if vC:
+                    if iE < bufEvent:
+                        bufX.append(vX)
+                        bufY.append(vY)                        
+                    else:
+                        h.Fill(vX,vY)
+            except:
+                pass
+                            
+            # if not a collection escape from loop
+            if not isCollection:
+                break
+
+        # create Histogram
+        if (iE+1) == bufEvent:
+            if len(bufX)==0:
+                minX=0
+            else:    
+                minX = min(bufX)
+                
+            if minX<0:
+                minX *= 1.1
+            elif minX>0:
+                minX *= 0.9
+            else:
+                minX = -1
+                
+            if len(bufX)==0:
+                maxX=0
+            else:    
+                maxX = max(bufX)
+                
+            if maxX<0:
+                maxX *= 0.9
+            elif maxX>0:
+                maxX *= 1.1
+            else:
+                maxX = 1
+
+            # create histogram if hist is None
+            if hist == None:
+                lpath = '/stat/tmp/PyKernelHist'
+                unregister(lpath)
+                # determine title of histo
+                if callable(valueX):
+                    titleX = valueX.__name__
+                else:
+                    titleX = valueX                    
+                if callable(valueY):
+                    titleY = valueY.__name__
+                else:
+                    titleY = valueY                    
+                h = bookProf(lpath, titleY+" vs "+titleX, 100, minX, maxX)
+            else:
+                h = hist
+
+    # buffered elements
+    for iB in range(len(bufX)):
+        vX = bufX[iB]
+        vY = bufY[iB]
+        h.Fill(vX,vY)
+        
+    return h    
+
+
+# plot profileD histogram
+def plotProf (classAndKey, valueX="$x", valueY="$x", criteria="True", nEvent=100):
+    """Plot profile-histogram
+    
+    :param classAndKey: combination of class name and key separeted with '#', 'Class#Key'
+    :param valueX: physics parameter for X in string
+    :param valueY: physics parameter for Y in string
+    :param criteria: selection criteria
+    :param nEvent: number of event to be processed
+
+    For detail, see `plot`
+    
+    """
+    h = fillProf (None, classAndKey, valueX, valueY, criteria, nEvent)
+    # draw 
+    h.Draw()
+    # this return is needed to draw up the canvas.
+    # note that PyKernel._hSave = h doesn't work although it keeps the hist in memory  
+    return h
+
+    
+# parse string
+def _parseString (str):
+    # remove $
+    str = re.sub("\$", "", str)
+    # replace XXX#YYY with StoreGate access
+    cK = re.findall(r'([\w_]+)#([\w_\*]+)',str)
+    for iCK in cK:
+        # when XXX#*
+        if iCK[1]=='*':
+            bStr = iCK[0]+"#\*"
+            aStr = 'retrieve(GNS.'+iCK[0]+')'
+        else:
+            bStr = iCK[0]+"#"+iCK[1]
+            aStr = 'retrieve(GNS.'+iCK[0]+',"'+iCK[1]+'")'
+        str = re.sub(bStr, aStr, str)
+    return str    
+
+
+# dump objects in StoreGate
+def dumpSG ():
+    '''
+    Dump objects in StoreGate
+
+    **examples**::
+
+      athena> dumpSG()
+      
+    '''
+    print GNS.StoreGate.pointer().dump()
+
+
+# book 1D/2D-histogram
+def book (*args):
+    '''
+    Book AIDA histogram (1D/2D) and register it to HistogramSvc.
+    Return reference to the histogram
+    
+    for 1D-histogram
+      book (path, title, nbinx, xmin, xmax)
+             
+    for 2D-histogram
+      book (path, title, nbinx, xmin, xmax, nbiny, ymin, ymax)
+
+    :Parameters:
+      - `args`: variable length parameter
+      
+             * path : path to the histogram
+
+             * title : title property of the histogram
+             
+             * nbinx : number of bins on the axis X
+             
+             * xmin : lower histogram edge on the axis X
+              
+             * xmax : upper histogram edge on the axis X
+
+             * nbiny : number of bins on the axis Y
+             
+             * ymin : lower histogram edge on the axis Y
+             
+             * ymax : upper histogram edge on the axis Y
+
+    **examples**::
+
+      athena> h = book("/stat/test1","Test1",100,0,100*GeV)
+      athena> h = book("/stat/test2","Test2",100,0,100*GeV,100,-50*GeV,50*GeV)
+
+    For detail, see `PyAnalysisExamples/HistTest.py`_
+
+    .. _PyAnalysisExamples/HistTest.py: http://atlas-sw.cern.ch/cgi-bin/viewcvs-atlas.cgi/offline/PhysicsAnalysis/PyAnalysis/PyAnalysisExamples/share/HistTest.py?rev=HEAD&content-type=text/vnd.viewcvs-markup
+      
+    '''
+    h = apply(theApp.histSvc().book, args)
+    # 1D hist
+    if (str(h).find('Histogram 1D') != -1):
+        return PyKHist.PyKHist1D(h)
+    # 2D hist
+    if (str(h).find('Histogram 2D') != -1):
+        return PyKHist.PyKHist2D(h)
+    # else
+    return None
+
+
+# book Profile-histogram
+def bookProf (*args):
+    '''
+    Book AIDA profile-histogram and register it to HistogramSvc.
+    Return reference to the histogram
+
+    Usually
+      bookProf (path, title, nbinx, xmin, xmax)
+
+    :Parameters:
+      - `args`: variable length parameter
+      
+             * path : path to the histogram
+
+             * title : title property of the histogram
+             
+             * nbinx : number of bins on the axis X
+             
+             * xmin : lower histogram edge on the axis X
+              
+             * xmax : upper histogram edge on the axis X
+
+    **examples**::
+
+      athena> h = bookProf("/stat/testp","Testp",100,0,100*GeV)
+
+    For detail, see `book`
+      
+    '''
+    return PyKHist.PyKHistProf(apply(theApp.histSvc()._ihs.bookProf, args))
+
+
+# unregister histogram from HistogramSvc
+def unregister (path):
+    '''
+    Unregister histogram from HistogramSvc
+    
+    :param path: path to the histogram
+
+    **examples**::
+
+      athena> unregister("/stat/tmpHist")
+      
+    '''
+    return theApp.histSvc().unregisterObject(path)
+
+
+# dump histograms in HistogramSvc
+def dumpHist ():
+    '''
+    Dump histograms in HistogramSvc
+
+    **examples**::
+
+      athena> dumpHist()
+      
+    '''
+    theApp.histSvc().dump()
+
+
+# retrieve histograms from HistogramSvc
+def retrieveHist (path):
+    '''
+    Retrieve histograms from HistogramSvc
+
+    :param path: path to the histogram
+
+    **examples**::
+
+      athena> h = retrieveHist("/stat/elec_n")
+      
+    '''
+    h = theApp.histSvc().retrieve(path)
+    # 1D hist
+    if (str(h).find('Histogram 1D') != -1):
+        return PyKHist.PyKHist1D(h)
+    # 2D hist
+    if (str(h).find('Histogram 2D') != -1):
+        return PyKHist.PyKHist2D(h)
+    # else
+    return None
+
+
+# set event loop type to pre-process
+def preProcess ():
+    '''
+    Set event loop type to pre-process
+
+    '''
+    global curEvent
+    curEvent = theApp.curEvent()
+    global eventLoopType    
+    eventLoopType = _PreProcess
+
+
+# set event loop type to normal-process
+def normalProcess ():
+    '''
+    Set event loop type to normal-process
+
+    '''
+    global eventLoopType    
+    eventLoopType = _NormalProcess
+    
+
+# set event loop type to hybrid-process
+def hybridProcess ():
+    '''
+    Set event loop type to hybrid-process
+
+    '''
+    global eventLoopType
+    eventLoopType = _HybridProcess
diff --git a/Control/PyKernel/python/__init__.py b/Control/PyKernel/python/__init__.py
new file mode 100755
index 0000000000000000000000000000000000000000..74583d364ec2ca794156596c7254d9b234a940c6
--- /dev/null
+++ b/Control/PyKernel/python/__init__.py
@@ -0,0 +1,2 @@
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
diff --git a/Control/PyKernel/share/InitPyKernel.py b/Control/PyKernel/share/InitPyKernel.py
new file mode 100755
index 0000000000000000000000000000000000000000..48af813f5c30cac79f9470f805bd91c11b4fe166
--- /dev/null
+++ b/Control/PyKernel/share/InitPyKernel.py
@@ -0,0 +1,30 @@
+from PyKernel import PyKernel
+
+from PyKernel.PyKernel import plot,plot2,fill,fill2,dumpSG,book,bookProf
+from PyKernel.PyKernel import plotProf,fillProf,unregister,dumpHist,retrieveHist
+
+import PyCintex
+PyCintex.loadDictionary("libPyKernelDict")
+
+# Bind the C++ global namespace to the Python variable g 
+g = PyCintex.Namespace('')
+
+## temporary hack to 'fix' #58078
+if 0:
+    # set ROOT stream
+    import re
+    import ROOT
+    rootStream=None
+    from AthenaCommon.AppMgr import ServiceMgr as svcMgr
+    EventSelector = svcMgr.EventSelector
+    if hasattr(EventSelector,"CollectionType") and  EventSelector.CollectionType == "ExplicitROOT":
+        filename = EventSelector.InputCollections[0]
+        if not re.search('\.root$',filename):
+            filename = filename+'.root'
+        file = ROOT.TFile(filename)
+        rootStream=ROOT.gROOT.FindObject('CollectionTree')
+
+    # this must be the last one
+    PyKernel.init(theApp,rootStream)
+
+    pass # temporary hack
diff --git a/Control/PyKernel/src/PyReverseProxy.cxx b/Control/PyKernel/src/PyReverseProxy.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..a28b51682a94428d6e58091fad4b8c3a62760340
--- /dev/null
+++ b/Control/PyKernel/src/PyReverseProxy.cxx
@@ -0,0 +1,8 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "PyKernel/PyReverseProxy.h"
+
+std::map<std::string, PyReverseProxy *> PyReverseProxy::m_proxyMap
+= std::map<std::string, PyReverseProxy *>();