diff --git a/pymo/preprocessing.py b/pymo/preprocessing.py
index 01f3a54..a7e12e6 100644
--- a/pymo/preprocessing.py
+++ b/pymo/preprocessing.py
@@ -96,6 +96,10 @@ def _to_pos(self, X):
                                      f[1]['%s_Yrotation'%joint], 
                                      f[1]['%s_Zrotation'%joint]] for f in rc.iterrows()]
 
+                    ################# in euler angle, the order of rotation axis is very important #####################
+                    rotation_order = rc.columns[0][rc.columns[0].find('rotation') - 1] + rc.columns[1][rc.columns[1].find('rotation') - 1] + rc.columns[2][rc.columns[2].find('rotation') - 1] #rotation_order is string : 'XYZ' or'ZYX' or ...
+                    ####################################################################################################
+
                 if pc.shape[1] < 3:
                     pos_values = [[0,0,0] for f in pc.iterrows()]
                 else:
@@ -107,7 +111,9 @@ def _to_pos(self, X):
                 #pos_values = [[0,0,0] for f in pc.iterrows()] #for deugging
                 
                 # Convert the eulers to rotation matrices
-                rotmats = np.asarray([Rotation([f[0], f[1], f[2]], 'euler', from_deg=True).rotmat for f in euler_values])                  
+                ############################ input rotation order as Rotation class's argument #########################
+                rotmats = np.asarray([Rotation([f[0], f[1], f[2]], 'euler', rotation_order, from_deg=True).rotmat for f in euler_values])
+                ########################################################################################################
                 tree_data[joint]=[
                                     [], # to store the rotation matrix
                                     []  # to store the calculated position
diff --git a/pymo/rotation_tools.py b/pymo/rotation_tools.py
index 1fe6072..ba208a9 100644
--- a/pymo/rotation_tools.py
+++ b/pymo/rotation_tools.py
@@ -18,8 +18,9 @@ def rad2deg(x):
     return x/math.pi*180
 
 class Rotation():
-    def __init__(self,rot, param_type, **params):
+    def __init__(self,rot, param_type, rotation_order, **params):
         self.rotmat = []
+        self.rotation_order = rotation_order
         if param_type == 'euler':
             self._from_euler(rot[0],rot[1],rot[2], params)
         elif param_type == 'expmap':
@@ -55,10 +56,15 @@ def _from_euler(self, alpha, beta, gamma, params):
 
         self.rotmat = np.eye(3)
 
-        self.rotmat = np.matmul(Rx, self.rotmat)
-        self.rotmat = np.matmul(Ry, self.rotmat)
-        self.rotmat = np.matmul(Rz, self.rotmat)
-        # self.rotmat = np.matmul(np.matmul(Rz, Ry), Rx)
+        ############################ inner product rotation matrix in order defined at BVH file #########################
+        for axis in self.rotation_order :
+            if axis == 'X' :
+                self.rotmat = np.matmul(Rx, self.rotmat)
+            elif axis == 'Y':
+                self.rotmat = np.matmul(Ry, self.rotmat)
+            else :
+                self.rotmat = np.matmul(Rz, self.rotmat)
+        ################################################################################################################
    
     def _from_expmap(self, alpha, beta, gamma, params):
         if (alpha == 0 and beta == 0 and gamma == 0):