diff --git a/src/block.cc b/src/block.cc
index ba255b7..930205d 100644
--- a/src/block.cc
+++ b/src/block.cc
@@ -120,6 +120,12 @@ namespace MultiArrayHelper
 	return *this;
     }
 
+    template <typename T>
+    size_t Block<T>::stepSize() const
+    {
+	return 1;
+    }
+    
     /**************
      *   MBlock   *
      **************/
@@ -157,6 +163,12 @@ namespace MultiArrayHelper
 	mBegPtr = &(*mData)[npos];
 	return *this;
     }   
+
+    template <typename T>
+    size_t MBlock<T>::stepSize() const
+    {
+	return 1;
+    }
     
     /******************
      *   BlockValue   *
@@ -167,7 +179,7 @@ namespace MultiArrayHelper
 			      size_t pos, size_t size) :
 	BlockBase<T>(size),
 	mData(&data),
-	mVal(data[pos]) {}
+	mVal(&data[pos]) {}
 
     template <typename T>
     BlockType BlockValue<T>::type() const
@@ -179,16 +191,22 @@ namespace MultiArrayHelper
     const T& BlockValue<T>::operator[](size_t i) const
     {
 	
-	return mVal;
+	return *mVal;
     }
 
     template <typename T>
     BlockValue<T>& BlockValue<T>::set(size_t npos)
     {
-	mVal = (*mData)[npos];
+	mVal = &(*mData)[npos];
 	return *this;
     }
 
+    template <typename T>
+    size_t BlockValue<T>::stepSize() const
+    {
+	return 0;
+    }
+    
     /*******************
      *   MBlockValue   *
      *******************/
@@ -198,7 +216,7 @@ namespace MultiArrayHelper
 				size_t pos, size_t size) :
 	MutableBlockBase<T>(size),
 	mData(&data),
-	mVal(data[pos]) {}
+	mVal(&data[pos]) {}
 
     template <typename T>
     BlockType MBlockValue<T>::type() const
@@ -210,23 +228,29 @@ namespace MultiArrayHelper
     const T& MBlockValue<T>::operator[](size_t i) const
     {
 	
-	return mVal;
+	return *mVal;
     }
 
     template <typename T>
     T& MBlockValue<T>::operator[](size_t i)
     {
 	
-	return mVal;
+	return *mVal;
     }
     
     template <typename T>
     MBlockValue<T>& MBlockValue<T>::set(size_t npos)
     {
-	mVal = (*mData)[npos];
+	mVal = &(*mData)[npos];
 	return *this;
     }
 
+    template <typename T>
+    size_t MBlockValue<T>::stepSize() const
+    {
+	return 0;
+    }
+    
     /******************
      *   SplitBlock   *
      ******************/
@@ -235,9 +259,9 @@ namespace MultiArrayHelper
     SplitBlock<T>::SplitBlock(const std::vector<T>& data, size_t begPos,
 			      size_t stepSize, size_t size) :
 	BlockBase<T>(size),
-	mData(data),
+	mData(&data),
 	mStepSize(stepSize),
-	mBegPtr(data.data() + begPos) {}
+	mBegPtr(&data[begPos]) {}
 
     template <typename T>
     BlockType SplitBlock<T>::type() const
@@ -259,6 +283,12 @@ namespace MultiArrayHelper
 	return *this;
     }
 
+    template <typename T>
+    size_t SplitBlock<T>::stepSize() const
+    {
+	return mStepSize;
+    }
+    
     /*******************
      *   MSplitBlock   *
      *******************/
@@ -267,8 +297,9 @@ namespace MultiArrayHelper
     MSplitBlock<T>::MSplitBlock(std::vector<T>& data, size_t begPos,
 				size_t stepSize, size_t size) :
 	MutableBlockBase<T>(size),
+	mData(&data),
 	mStepSize(stepSize),
-	mBegPtr(data.data() + begPos) {}
+	mBegPtr(&data[begPos]) {}
 
     template <typename T>
     BlockType MSplitBlock<T>::type() const
@@ -296,6 +327,12 @@ namespace MultiArrayHelper
 	mBegPtr = &(*mData)[npos];
 	return *this;
     }
+
+    template <typename T>
+    size_t MSplitBlock<T>::stepSize() const
+    {
+	return mStepSize;
+    }
     
     /*******************
      *   BlockResult   *
@@ -331,5 +368,11 @@ namespace MultiArrayHelper
     {
 	return *this;
     }
+
+    template <typename T>
+    size_t BlockResult<T>::stepSize() const
+    {
+	return 1;
+    }
     
 } // end namespace MultiArrayHelper
diff --git a/src/block.h b/src/block.h
index c129402..ee6f7a7 100644
--- a/src/block.h
+++ b/src/block.h
@@ -40,6 +40,7 @@ namespace MultiArrayHelper
 	BlockBase& operator=(BlockBase&& res) = default;
 	
 	virtual BlockType type() const = 0;
+	virtual size_t stepSize() const = 0;
 	
 	virtual size_t size() const;
 	virtual const T& operator[](size_t pos) const = 0;
@@ -58,6 +59,16 @@ namespace MultiArrayHelper
 	size_t mSize;
     };
 
+    template <typename T>
+    std::ostream& operator<<(std::ostream& out, const BlockBase<T>& block)
+    {
+	out << block[0];
+	for(size_t i = 1; i != block.size(); ++i){
+	    out << ", " << block[i];
+	}
+	return out;
+    }
+    
     template <typename T>
     class MutableBlockBase : public BlockBase<T>
     {
@@ -85,6 +96,7 @@ namespace MultiArrayHelper
 	virtual BlockType type() const override;
 	virtual const T& operator[](size_t pos) const override;
 	virtual Block& set(size_t npos) override;
+	virtual size_t stepSize() const override;
 	
     protected:
 	const std::vector<T>* mData;
@@ -102,6 +114,7 @@ namespace MultiArrayHelper
 	virtual const T& operator[](size_t pos) const override;
 	virtual T& operator[](size_t pos) override;
 	virtual MBlock& set(size_t npos) override;
+	virtual size_t stepSize() const override;
 	
     protected:
 	std::vector<T>* mData;
@@ -119,10 +132,11 @@ namespace MultiArrayHelper
 	virtual BlockType type() const override;
 	virtual const T& operator[](size_t pos) const override;
 	virtual BlockValue& set(size_t npos) override;
+	virtual size_t stepSize() const override;
 	
     protected:
 	const std::vector<T>* mData;
-	T& mVal;
+	const T* mVal;
     };
 
     template <typename T>
@@ -137,10 +151,11 @@ namespace MultiArrayHelper
 	virtual const T& operator[](size_t pos) const override;
 	virtual T& operator[](size_t pos) override;
 	virtual MBlockValue& set(size_t npos) override;
+	virtual size_t stepSize() const override;
 	
     protected:
 	std::vector<T>* mData;
-	T& mVal;
+	T* mVal;
     };
     
     template <typename T>
@@ -155,6 +170,7 @@ namespace MultiArrayHelper
 	virtual BlockType type() const override;
 	virtual const T& operator[](size_t pos) const override;
 	virtual SplitBlock& set(size_t npos) override;
+	virtual size_t stepSize() const override;
 	
     protected:
 	const std::vector<T>* mData;
@@ -175,6 +191,7 @@ namespace MultiArrayHelper
 	virtual const T& operator[](size_t pos) const override;
 	virtual T& operator[](size_t pos) override;
 	virtual MSplitBlock& set(size_t npos) override;
+	virtual size_t stepSize() const override;
 	
     protected:
 	std::vector<T>* mData;
@@ -197,6 +214,7 @@ namespace MultiArrayHelper
 	virtual const T& operator[](size_t pos) const override;
 	virtual T& operator[](size_t i) override;
 	virtual BlockResult& set(size_t npos) override;
+	virtual size_t stepSize() const override;
 	
     protected:
 	std::vector<T> mRes;
diff --git a/src/multi_array.cc b/src/multi_array.cc
index 8da67ff..58df5ec 100644
--- a/src/multi_array.cc
+++ b/src/multi_array.cc
@@ -408,6 +408,14 @@ namespace MultiArrayTools
 	return OperationRoot<T,SRanges...>(*this, inds...);
     }
 
+    template <typename T, class... SRanges>
+    ConstOperationRoot<T,SRanges...>
+    MutableMultiArrayBase<T,SRanges...>::operator()(std::shared_ptr<typename SRanges::IndexType>&... inds) const
+    {
+	return ConstOperationRoot<T,SRanges...>(*this, inds...);
+    }
+
+    
     /*******************
      *  MultiArray     *	     
      *******************/
diff --git a/src/multi_array.h b/src/multi_array.h
index 973ff12..07d6564 100644
--- a/src/multi_array.h
+++ b/src/multi_array.h
@@ -96,7 +96,8 @@ namespace MultiArrayTools
 
 	virtual bool isConst() const;
 
-	ConstOperationRoot<T,SRanges...> operator()(std::shared_ptr<typename SRanges::IndexType>&... inds) const;
+	virtual ConstOperationRoot<T,SRanges...>
+	operator()(std::shared_ptr<typename SRanges::IndexType>&... inds) const;
 	
 	virtual bool isInit() const;
 
@@ -179,7 +180,9 @@ namespace MultiArrayTools
 
 	virtual bool isConst() const override;
 
-	OperationRoot<T,SRanges...> operator()(std::shared_ptr<typename SRanges::IndexType>&... inds);
+	virtual ConstOperationRoot<T,SRanges...>
+	operator()(std::shared_ptr<typename SRanges::IndexType>&... inds) const override;
+	virtual OperationRoot<T,SRanges...> operator()(std::shared_ptr<typename SRanges::IndexType>&... inds);
     };
     
     template <typename T, class... SRanges>
diff --git a/src/multi_array_header.h b/src/multi_array_header.h
index 5fa3f05..6e37477 100644
--- a/src/multi_array_header.h
+++ b/src/multi_array_header.h
@@ -53,6 +53,13 @@ namespace MultiArrayTools
     typedef SingleRange<std::array<int,4>, RangeType::MOMENTUM> LinMom4dR;
     // ...
     */
+
+    template <typename... T>
+    std::ostream& operator<<(std::ostream& out, const std::tuple<T...>& tp)
+    {
+	PackNum<sizeof...(T)-1>::printTuple(out, tp);
+	return out;
+    }
 }
 
 #endif
diff --git a/src/multi_array_operation.cc b/src/multi_array_operation.cc
index d8f8525..a6b5cc1 100644
--- a/src/multi_array_operation.cc
+++ b/src/multi_array_operation.cc
@@ -28,10 +28,11 @@ namespace MultiArrayTools
 	// stepSize == 1 => BLOCK
 	// stepSize > 1 => SPLIT :)
 	BTSS out(BlockType::VALUE, 0);
-	for(size_t inum = 0; inum != i->rangePtr()->dim(); ++inum){
+	size_t lastNum = i->rangePtr()->dim();
+	for(size_t inum = 0; inum != lastNum; ++inum){
 	    auto ii = i->getPtr(inum);
 	    if(ii == j){
-		if(inum == 0 and first){
+		if(inum == lastNum - 1 and first){
 		    out = BTSS(BlockType::BLOCK, 1);
 		}
 		else {
@@ -42,7 +43,7 @@ namespace MultiArrayTools
 	    
 	    if(ii->type() == IndexType::MULTI or
 	       ii->type() == IndexType::CONT){	
-		BTSS tmp = getBlockType(ii, j, inum == 0);
+		BTSS tmp = getBlockType(ii, j, inum == lastNum - 1);
 		if(tmp.first != BlockType::VALUE){
 		    out = tmp;
 		}
@@ -54,7 +55,6 @@ namespace MultiArrayTools
     template <typename T>
     std::shared_ptr<BlockBase<T> > makeBlock(const std::vector<T>& vec, size_t stepSize, size_t blockSize)
     {
-	VCHECK(stepSize);
 	if(stepSize == 0){
 	    return std::make_shared<BlockValue<T> >(vec, 0, blockSize);
 	}
@@ -69,7 +69,6 @@ namespace MultiArrayTools
     template <typename T>
     std::shared_ptr<MutableBlockBase<T> > makeBlock(std::vector<T>& vec, size_t stepSize, size_t blockSize)
     {
-	VCHECK(stepSize);
 	if(stepSize == 0){
 	    return std::make_shared<MBlockValue<T> >(vec, 0, blockSize);
 	}
@@ -202,6 +201,7 @@ namespace MultiArrayTools
 	auto blockIndex = seekBlockIndex( mIndex, second);
 	block(blockIndex);
 	second.block(blockIndex);
+
 	for(*mIndex = 0; mIndex->pos() != mIndex->max(); mIndex->pp(blockIndex) ){
 	    get() = mSecond.get();
 	}
@@ -232,7 +232,7 @@ namespace MultiArrayTools
     template <typename T, class... Ranges>
     const OperationMaster<T,Ranges...>& OperationMaster<T,Ranges...>::block() const
     {
-	mBlockPtr->set( mIndex->pos() ); // containerRange <-> multiRange !!!
+	mBlockPtr->set( mIndex->pos() );
 	return *this;
     }
     
@@ -247,7 +247,7 @@ namespace MultiArrayTools
 	OperationBase<T>(), OperationTemplate<T,ConstOperationRoot<T,Ranges...> >(this),
 	mArrayRef(ma), mIndex( std::make_shared<IndexType>( mArrayRef.range() ) )
     {
-	mIndex(indices...);
+	(*mIndex)(indices...);
     }
 
     template <typename T, class... Ranges>
@@ -268,14 +268,14 @@ namespace MultiArrayTools
     template <typename T, class... Ranges>
     const ConstOperationRoot<T,Ranges...>& ConstOperationRoot<T,Ranges...>::block() const
     {
-	mBlockPtr->set( mIndex->pos() );
+	mBlockPtr->set( (*mIndex)().pos() );
 	return *this;
     }
     
     /***********************
      *   OperationRoot     *
      ***********************/
-
+    
     template <typename T, class... Ranges>
     OperationRoot<T,Ranges...>::
     OperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
@@ -298,7 +298,7 @@ namespace MultiArrayTools
 	block();
 	return *mBlockPtr;
     }
-
+    
     template <typename T, class... Ranges>
     MutableBlockBase<T>& OperationRoot<T,Ranges...>::get()
     {
@@ -317,13 +317,13 @@ namespace MultiArrayTools
     template <typename T, class... Ranges>
     const OperationRoot<T,Ranges...>& OperationRoot<T,Ranges...>::block() const
     {
-	mBlockPtr->set( mIndex->pos() );
+	mBlockPtr->set( (*mIndex)().pos() );
 	return *this;
     }
     
-    /***********************
-     *   OperationRoot     *
-     ***********************/
+    /*******************
+     *   Operation     *
+     *******************/
     
     template <typename T, class OpFunction, class... Ops>
     Operation<T,OpFunction,Ops...>::Operation(const Ops&... ops) :
diff --git a/src/multi_array_operation.h b/src/multi_array_operation.h
index 276a755..8e2f6a0 100644
--- a/src/multi_array_operation.h
+++ b/src/multi_array_operation.h
@@ -167,7 +167,7 @@ namespace MultiArrayTools
 	std::shared_ptr<IndexType> mIndex;
  	mutable std::shared_ptr<BlockBase<T> > mBlockPtr;
     };
-
+    
     template <typename T, class... Ranges>
     class OperationRoot : public MutableOperationBase<T>,
 			  public OperationTemplate<T,OperationRoot<T,Ranges...> >
@@ -197,7 +197,7 @@ namespace MultiArrayTools
 	std::shared_ptr<IndexType> mIndex;
 	mutable std::shared_ptr<MutableBlockBase<T> > mBlockPtr;
     };
-
+    
     template <typename T, class OpFunction, class... Ops>
     class Operation : public OperationBase<T>,
 		      public OperationTemplate<T,Operation<T,OpFunction,Ops...> >
diff --git a/src/op_unit_test.cc b/src/op_unit_test.cc
index d743049..804347c 100644
--- a/src/op_unit_test.cc
+++ b/src/op_unit_test.cc
@@ -114,8 +114,8 @@ namespace {
     TEST_F(OpTest_MDim, ExecOp1)
     {
 	MultiArray<double,SRange,SRange> res(sr2ptr,sr4ptr);
-	MultiArray<double,SRange> ma1(sr2ptr, v1);
-	MultiArray<double,SRange> ma2(sr4ptr, v2);
+	const MultiArray<double,SRange> ma1(sr2ptr, v1);
+	const MultiArray<double,SRange> ma2(sr4ptr, v2);
 
 	auto i1 = std::dynamic_pointer_cast<SRange::IndexType>( sr2ptr->index() );
 	auto i2 = std::dynamic_pointer_cast<SRange::IndexType>( sr4ptr->index() );
diff --git a/src/pack_num.h b/src/pack_num.h
index ba91b6b..c5279b2 100644
--- a/src/pack_num.h
+++ b/src/pack_num.h
@@ -231,6 +231,12 @@ namespace MultiArrayHelper
 	    btv.insert(btv.end(), subvec.begin(), subvec.end() );
 	    PackNum<N-1>::makeBlockTypeVec(btv, ops, idxPtr);
 	}
+
+	template <typename... T>
+	static void printTuple(std::ostream& out, const std::tuple<T...>& tp){
+	    out << std::get<sizeof...(T)-N-1>(tp) << ", ";
+	    PackNum<N-1>::printTuple(out, tp);
+	}
 	
     };
     
@@ -403,7 +409,12 @@ namespace MultiArrayHelper
 	    auto subvec = std::move( std::get<0>(ops).block(idxPtr) );
 	    btv.insert(btv.end(), subvec.begin(), subvec.end() );
 	}
-	
+
+	template <typename... T>
+	static void printTuple(std::ostream& out, const std::tuple<T...>& tp){
+	    out << std::get<sizeof...(T)-1>(tp);
+	}
+
     };
     
 } // end namespace MultiArrayHelper