diff --git a/programs/Hadrons/Application.cc b/programs/Hadrons/Application.cc
index 923bd0eb..8f132056 100644
--- a/programs/Hadrons/Application.cc
+++ b/programs/Hadrons/Application.cc
@@ -59,6 +59,8 @@ void Application::run(void)
     
     g.addEdge("A", "B");
     g.addEdge("B", "D");
+    g.addEdge("B", "X1");
+    g.addEdge("X1", "X2");
     g.addEdge("D", "E");
     g.addEdge("E", "C");
     g.addEdge("Z", "Y");
@@ -75,13 +77,16 @@ void Application::run(void)
     }
     for (auto &h: vec)
     {
-        auto top = h.topoSort();
-        while (!top.empty())
+        auto top = h.allTopoSort();
+        for (auto &s: top)
         {
-            cout << top.top() << " ";
-            top.pop();
+            for (auto &v: s)
+            {
+                cout << v << " ";
+            }
+            cout << endl;
         }
-        cout << endl;
+        cout << "--------" << endl;
     }
 }
 
diff --git a/programs/Hadrons/Graph.hpp b/programs/Hadrons/Graph.hpp
index 15eba21a..14077ed0 100644
--- a/programs/Hadrons/Graph.hpp
+++ b/programs/Hadrons/Graph.hpp
@@ -58,12 +58,13 @@ public:
     // tests
     bool gotValue(const T &value) const;
     // graph topological manipulations
-    std::vector<T>        getAdjacentVertices(const T &value) const;
-    std::vector<T>        getChildren(const T &value) const;
-    std::vector<T>        getParents(const T &value) const;
-    std::vector<T>        getRoots(void) const;
-    std::vector<Graph<T>> getConnectedComponents(void) const;
-    std::stack<T>         topoSort(void);
+    std::vector<T>              getAdjacentVertices(const T &value) const;
+    std::vector<T>              getChildren(const T &value) const;
+    std::vector<T>              getParents(const T &value) const;
+    std::vector<T>              getRoots(void) const;
+    std::vector<Graph<T>>       getConnectedComponents(void) const;
+    std::vector<T>              topoSort(void);
+    std::vector<std::vector<T>> allTopoSort(void);
     // I/O
     friend std::ostream & operator<<(std::ostream &out, const Graph<T> &g)
     {
@@ -96,8 +97,8 @@ private:
     void depthFirstSearch(void);
     void depthFirstSearch(const T &root);
 private:
-    std::map<T, bool> isMarked_;
-    std::set<Edge>    edgeSet_;
+    std::map<T, bool>              isMarked_;
+    std::set<Edge>                 edgeSet_;
 };
 
 /******************************************************************************
@@ -422,12 +423,13 @@ std::vector<Graph<T>> Graph<T>::getConnectedComponents(void) const
 
 // topological sort using Tarjan's algorithm
 template <typename T>
-std::stack<T> Graph<T>::topoSort(void)
+std::vector<T> Graph<T>::topoSort(void)
 {
-    std::stack<T>     res;
+    std::stack<T>     buf;
+    std::vector<T>    res;
     const T           *vPt;
     std::map<T, bool> tmpMarked(isMarked_);
-    
+
     // visit function
     std::function<void(const T &)> visit = [&](const T &v)
     {
@@ -446,7 +448,7 @@ std::stack<T> Graph<T>::topoSort(void)
             }
             this->mark(v);
             tmpMarked[v] = false;
-            res.push(v);
+            buf.push(v);
         }
     };
     
@@ -455,6 +457,7 @@ std::stack<T> Graph<T>::topoSort(void)
     {
         tmpMarked.at(v.first) = false;
     }
+    
     // loop on unmarked vertices
     vPt = getFirstUnmarked();
     while (vPt)
@@ -464,6 +467,92 @@ std::stack<T> Graph<T>::topoSort(void)
     }
     unmarkAll();
     
+    // create result vector
+    while (!buf.empty())
+    {
+        res.push_back(buf.top());
+        buf.pop();
+    }
+    
+    return res;
+}
+
+// generate all possible topological sorts
+// Y. L. Varol & D. Rotem, Comput. J. 24(1), pp. 83–84, 1981
+// http://comjnl.oupjournals.org/cgi/doi/10.1093/comjnl/24.1.83
+template <typename T>
+std::vector<std::vector<T>> Graph<T>::allTopoSort(void)
+{
+    std::vector<std::vector<T>>    res;
+    std::map<T, std::map<T, bool>> iMat;
+    
+    // create incidence matrix
+    for (auto &v1: isMarked_)
+    for (auto &v2: isMarked_)
+    {
+        iMat[v1.first][v2.first] = false;
+    }
+    for (auto &v: isMarked_)
+    {
+        auto cVec = getChildren(v.first);
+        
+        for (auto &c: cVec)
+        {
+            iMat[v.first][c] = true;
+        }
+    }
+    
+    // generate initial topological sort
+    res.push_back(topoSort());
+    
+    // generate all other topological sorts by permutation
+    std::vector<T>            p = res[0];
+    const unsigned int        n = size();
+    std::vector<unsigned int> loc(n);
+    unsigned int              i, k, k1;
+    T                         obj_k, obj_k1;
+    bool                      isFinal;
+    
+    for (unsigned int j = 0; j < n; ++j)
+    {
+        loc[j] = j;
+    }
+    i = 0;
+    while (i < n-1)
+    {
+        k      = loc[i];
+        k1     = k + 1;
+        obj_k  = p[k];
+        if (k1 >= n)
+        {
+            isFinal = true;
+            obj_k1  = obj_k;
+        }
+        else
+        {
+            isFinal = false;
+            obj_k1  = p[k1];
+        }
+        if (iMat[res[0][i]][obj_k1] or isFinal)
+        {
+            for (unsigned int l = k; l >= i + 1; --l)
+            {
+                p[l]   = p[l-1];
+            }
+            p[i]   = obj_k;
+            loc[i] = i;
+            i++;
+        }
+        else
+        {
+            p[k]   = obj_k1;
+            p[k1]  = obj_k;
+            loc[i] = k1;
+            i      = 0;
+            res.push_back(p);
+        }
+    }
+    
     return res;
 }