diff --git a/src/base_def.h b/src/base_def.h
index 6409a6e..cd71a3c 100644
--- a/src/base_def.h
+++ b/src/base_def.h
@@ -114,7 +114,7 @@ namespace MultiArrayTools
     class MultiArrayOperation;
 
     // multi_array_operation.h
-    template <typename T, class Operation, class Range, class MAOp>
+    template <typename T, class Operation, class Range, class... MAOps>
     class MultiArrayContraction;
     
     // slice.h
diff --git a/src/index_base.h b/src/index_base.h
index 67b5e2a..c766b17 100644
--- a/src/index_base.h
+++ b/src/index_base.h
@@ -40,7 +40,7 @@ namespace MultiArrayTools
 	virtual MultiRangeType rangeType() const = 0;
 
 	virtual bool link(IndefinitIndexBase* toLink);  
-	virtual void freeLinked();
+	virtual void freeLinked(); 
 	virtual bool linked() const;
 	virtual bool linkedTo(IndefinitIndexBase* link) const;
 	virtual void linkTo(IndefinitIndexBase* target) = 0; 
diff --git a/src/multi_array_operation.cc b/src/multi_array_operation.cc
index 113605a..fb9c25a 100644
--- a/src/multi_array_operation.cc
+++ b/src/multi_array_operation.cc
@@ -22,7 +22,7 @@ namespace MultiArrayTools
     }
 
     template <typename T>
-    void MultiArrayOperationBase<T>::freeIndex()
+    void MultiArrayOperationBase<T>::freeIndex() const
     {
 	mIibPtr->freeLinked();
     }
@@ -168,27 +168,37 @@ namespace MultiArrayTools
     }
 
     template <typename T, class Range>
-    template <class Range2, class ContractOperation>
-    MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperationRoot<T,Range> >
+    template <class Range2, class ContractOperation, class... MAOps>
+    MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperationRoot<T,Range>,MAOps...>
     MultiArrayOperationRoot<T,Range>::contract(const ContractOperation& cop,
-					       const std::string& indexName) const
+					       const std::string& indexName,
+					       const MAOps&... mao) const
     {
+	//typename Range2::IndexType ind;
+	//ind.name(indexName);
+
 	typename Range2::IndexType* ind = dynamic_cast<typename Range2::IndexType*>( mIndex.getLinked(indexName) );
 	return MultiArrayContraction<T,ContractOperation,Range2,
-				     MultiArrayOperationRoot<T,Range> >(cop, *this, *ind);
+				     MultiArrayOperationRoot<T,Range>,
+				     MAOps...>(cop, *ind, *this, mao...);
     }
 
     template <typename T, class Range>
-    template <class Range2, class ContractOperation>
-    MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperationRoot<T,Range> >
+    template <class Range2, class ContractOperation, class... MAOps>
+    MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperationRoot<T,Range>,MAOps...>
     MultiArrayOperationRoot<T,Range>::contract(const ContractOperation& cop,
 					       const std::string& indexName,
 					       const typename Range2::IndexType& begin,
-					       const typename Range2::IndexType& end) const
+					       const typename Range2::IndexType& end,
+					       const MAOps&... mao) const
     {
+	//typename Range2::IndexType ind;
+	//ind.name(indexName);
+
 	typename Range2::IndexType* ind = dynamic_cast<typename Range2::IndexType*>( mIndex.getLinked(indexName) );
 	return MultiArrayContraction<T,ContractOperation,Range2,
-				     MultiArrayOperationRoot<T,Range> >(cop, *this, *ind, begin, end);
+				     MultiArrayOperationRoot<T,Range>,
+				     MAOps...>(cop, *ind, begin, end, *this, mao...);
     }
 
     
@@ -324,7 +334,7 @@ namespace MultiArrayTools
     }
 
     template <typename T, class Range>
-    void MultiArrayOperationRoot<T,Range>::freeIndex()
+    void MultiArrayOperationRoot<T,Range>::freeIndex() const
     {
 	MAOB::mIibPtr->freeLinked();
 	mIndex = mArrayRef.beginIndex();
@@ -438,27 +448,31 @@ namespace MultiArrayTools
     }
 
     template <typename T, class Range>
-    template <class Range2, class ContractOperation>
-    MultiArrayContraction<T,ContractOperation,Range2,ConstMultiArrayOperationRoot<T,Range> >
+    template <class Range2, class ContractOperation, class... MAOps>
+    MultiArrayContraction<T,ContractOperation,Range2,ConstMultiArrayOperationRoot<T,Range>,MAOps...>
     ConstMultiArrayOperationRoot<T,Range>::contract(const ContractOperation& cop,
-						    const std::string& indexName) const
+						    const std::string& indexName,
+						    const MAOps&... mao) const
     {
 	typename Range2::IndexType* ind = dynamic_cast<typename Range2::IndexType*>( mIndex.getLinked(indexName) );
 	return MultiArrayContraction<T,ContractOperation,Range2,
-				     ConstMultiArrayOperationRoot<T,Range> >(cop, *this, *ind);
+				     ConstMultiArrayOperationRoot<T,Range>,
+				     MAOps...>(cop, *ind, *this, mao...);
     }
 
     template <typename T, class Range>
-    template <class Range2, class ContractOperation>
-    MultiArrayContraction<T,ContractOperation,Range2,ConstMultiArrayOperationRoot<T,Range> >
+    template <class Range2, class ContractOperation, class... MAOps>
+    MultiArrayContraction<T,ContractOperation,Range2,ConstMultiArrayOperationRoot<T,Range>,MAOps...>
     ConstMultiArrayOperationRoot<T,Range>::contract(const ContractOperation& cop,
 						    const std::string& indexName,
 						    const typename Range2::IndexType& begin,
-						    const typename Range2::IndexType& end) const
+						    const typename Range2::IndexType& end,
+						    const MAOps&... mao) const
     {
 	typename Range2::IndexType* ind = dynamic_cast<typename Range2::IndexType*>( mIndex.getLinked(indexName) );
 	return MultiArrayContraction<T,ContractOperation,Range2,
-				     ConstMultiArrayOperationRoot<T,Range> >(cop, *this, *ind, begin, end);
+				     ConstMultiArrayOperationRoot<T,Range>,
+				     MAOps...>(cop, *ind, begin, end, *this, mao...);
     }
 
 
@@ -546,7 +560,7 @@ namespace MultiArrayTools
     }
 
     template <typename T, class Range>
-    void ConstMultiArrayOperationRoot<T,Range>::freeIndex()
+    void ConstMultiArrayOperationRoot<T,Range>::freeIndex() const
     {
 	MAOB::mIibPtr->freeLinked();
 	mIndex = mArrayRef.beginIndex();
@@ -646,6 +660,38 @@ namespace MultiArrayTools
 				   MAOps2...>(op, *this, secs...);
     }
 
+    /*
+    template <typename T, class Operation, class... MAOps>
+    template <class Range2, class ContractOperation>
+    MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperation<T,Operation,MAOps...> >
+    MultiArrayOperation<T,Operation,MAOps...>::
+    contract(const ContractOperation& cop, const std::string& indexName) const
+    {
+	typename Range2::IndexType* ind = dynamic_cast<typename Range2::IndexType*>( mIndex.getLinked(indexName) );
+	//typename Range2::IndexType ind = Range2().begin();
+	//ind.name(indexName);
+	return MultiArrayContraction<T,ContractOperation,Range2,
+				     MultiArrayOperation<T,Operation,MAOps...> >(cop, *this, *ind);
+	
+    }
+    
+    template <typename T, class Operation, class... MAOps>
+    template <class Range2, class ContractOperation>
+    MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperation<T,Operation,MAOps...> >
+    MultiArrayOperation<T,Operation,MAOps...>::
+    contract(const ContractOperation& cop, const std::string& indexName,
+	     const typename Range2::IndexType& begin,
+	     const typename Range2::IndexType& end) const
+    {
+	typename Range2::IndexType* ind = dynamic_cast<typename Range2::IndexType*>( mIndex.getLinked(indexName) );
+	//typename Range2::IndexType ind = Range2().begin();
+	//ind.name(indexName);
+	return MultiArrayContraction<T,ContractOperation,Range2,
+				     MultiArrayOperation<T,Operation,MAOps...> >(cop, *this, *ind, begin, end);
+
+    }
+    */
+    
     template <typename T, class Operation, class... MAOps>
     template <class MAOp2>
     auto MultiArrayOperation<T,Operation,MAOps...>::operator+(const MAOp2& sec)
@@ -703,11 +749,12 @@ namespace MultiArrayTools
      *   MultiArrayContraction     *
      *******************************/
 
-    template <typename T, class ContractOperation, class Range, class MAOp>
-    MultiArrayContraction<T,ContractOperation,Range,MAOp>::
-    MultiArrayContraction(const ContractOperation& cop, const MAOp& mao,
-			  const typename Range::IndexType& runIndex) :
-	MultiArrayOperation<T,ContractOperation,MAOp>(cop, mao),
+    template <typename T, class ContractOperation, class Range, class... MAOps>
+    MultiArrayContraction<T,ContractOperation,Range,MAOps...>::
+    MultiArrayContraction(const ContractOperation& cop,
+			  const typename Range::IndexType& runIndex,
+			  const MAOps&... mao) :
+	MultiArrayOperation<T,ContractOperation,MAOps...>(cop, mao...),
 	mBeginIndex(runIndex), mEndIndex(runIndex),
 	mRunIndex(runIndex)
     {
@@ -716,13 +763,14 @@ namespace MultiArrayTools
 	MAO::linkIndicesTo(&mRunIndex);
     }
 
-    template <typename T, class ContractOperation, class Range, class MAOp>
-    MultiArrayContraction<T,ContractOperation,Range,MAOp>::
-    MultiArrayContraction(const ContractOperation& cop, const MAOp& mao,
+    template <typename T, class ContractOperation, class Range, class... MAOps>
+    MultiArrayContraction<T,ContractOperation,Range,MAOps...>::
+    MultiArrayContraction(const ContractOperation& cop,
 			  const typename Range::IndexType& runIndex,
 			  const typename Range::IndexType& beginIndex,
-			  const typename Range::IndexType& endIndex) :
-	MultiArrayOperation<T,ContractOperation,MAOp>(cop, mao),
+			  const typename Range::IndexType& endIndex,
+			  const MAOps&... mao) :
+	MultiArrayOperation<T,ContractOperation,MAOps...>(cop, mao...),
 	mBeginIndex(beginIndex), mEndIndex(endIndex),
 	mRunIndex(runIndex)
     {
@@ -731,14 +779,16 @@ namespace MultiArrayTools
 
     
     // for the moment simplest case only:
-    template <typename T, class ContractOperation, class Range, class MAOp>
-    const T& MultiArrayContraction<T,ContractOperation,Range,MAOp>::get() const
+    template <typename T, class ContractOperation, class Range, class... MAOps>
+    const T& MultiArrayContraction<T,ContractOperation,Range,MAOps...>::get() const
     {
 	MAO::mOp.reset();
 	for(mRunIndex.copyPos( mBeginIndex ); mRunIndex.pos() != mEndIndex.pos(); ++mRunIndex){
-	    MAO::mOp(std::get<0>(MAO::mArgs).get() );
+	    OperationCall<sizeof...(MAOps)-1>::
+		template callOperation(MAO::mOp, MAO::mArgs);
+	    //MAO::mOp(std::get<0>(MAO::mArgs).get() );
 	}
-	MAO::mOp.endOp(MAO::mVal);
+	MAO::mOp.endOp();
 	return MAO::mOp();
     }
 }
diff --git a/src/multi_array_operation.h b/src/multi_array_operation.h
index def9563..3cf184a 100644
--- a/src/multi_array_operation.h
+++ b/src/multi_array_operation.h
@@ -26,10 +26,10 @@ namespace MultiArrayTools
 	
 	virtual const T& get() const = 0;
 
-	virtual void freeIndex();
+	virtual void freeIndex() const;
 	
     protected:
-	IndefinitIndexBase* mIibPtr = nullptr;
+        mutable IndefinitIndexBase* mIibPtr = nullptr;
     };
     
     
@@ -74,15 +74,17 @@ namespace MultiArrayTools
 	MultiArrayOperation<T,Operation,MultiArrayOperationRoot<T,Range>, MAOps...>
 	operator()(const Operation& op, const MAOps&... secs) const;
 
-	template < class Range2, class ContractOperation>
-	MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperationRoot<T,Range> >
-	contract(const ContractOperation& cop, const std::string& indexName) const;
+	template < class Range2, class ContractOperation, class... MAOps>
+	MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperationRoot<T,Range>, MAOps...>
+	contract(const ContractOperation& cop, const std::string& indexName,
+		 const MAOps&... mao) const;
 	
-	template <class Range2, class ContractOperation>
-	MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperationRoot<T,Range> >
+	template <class Range2, class ContractOperation, class... MAOps>
+	MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperationRoot<T,Range>, MAOps...>
 	contract(const ContractOperation& cop, const std::string& indexName,
 		 const typename Range2::IndexType& begin,
-		 const typename Range2::IndexType& end) const;
+		 const typename Range2::IndexType& end,
+		 const MAOps&... mao) const;
 
 	
 	template <class MAOp>
@@ -128,7 +130,7 @@ namespace MultiArrayTools
 
 	const Name& name() const;
 
-	virtual void freeIndex() override;
+	virtual void freeIndex() const override;
 
 	const MultiArrayBase<T,Range>& getCont() const { return mArrayRef; }
 	
@@ -170,16 +172,18 @@ namespace MultiArrayTools
 	MultiArrayOperation<T,Operation,ConstMultiArrayOperationRoot<T,Range>, MAOps...>
 	operator()(const Operation& op, const MAOps&... secs) const;
 
-	template <class Range2, class ContractOperation>
-	MultiArrayContraction<T,ContractOperation,Range2,ConstMultiArrayOperationRoot<T,Range> >
-	contract(const ContractOperation& cop, const std::string& indexName) const;
+	template <class Range2, class ContractOperation, class... MAOps>
+	MultiArrayContraction<T,ContractOperation,Range2,ConstMultiArrayOperationRoot<T,Range>, MAOps...>
+	contract(const ContractOperation& cop, const std::string& indexName,
+		 const MAOps&... mao) const;
 
 	
-	template <class Range2, class ContractOperation>
-	MultiArrayContraction<T,ContractOperation,Range2,ConstMultiArrayOperationRoot<T,Range> >
+	template <class Range2, class ContractOperation, class... MAOps>
+	MultiArrayContraction<T,ContractOperation,Range2,ConstMultiArrayOperationRoot<T,Range>, MAOps...>
 	contract(const ContractOperation& cop, const std::string& indexName,
 		 const typename Range2::IndexType& begin,
-		 const typename Range2::IndexType& end) const;
+		 const typename Range2::IndexType& end,
+		 const MAOps&... mao) const;
 	
 	template <class MAOp>
 	auto operator+(const MAOp& sec) const -> decltype(operator()(std::plus<T>(), sec));
@@ -210,7 +214,7 @@ namespace MultiArrayTools
 
 	const Name& name() const;
 
-	virtual void freeIndex() override;
+	virtual void freeIndex() const override;
 
 	const MultiArrayBase<T,Range>& getCont() const { return mArrayRef; }
 	
@@ -253,6 +257,19 @@ namespace MultiArrayTools
 	template <class Operation2, class... MAOps2>
 	MultiArrayOperation<T,Operation2,MultiArrayOperation<T,Operation,MAOps...>,MAOps2...>
 	operator()(const Operation2& op, const MAOps2&... secs) const;
+
+	/*
+	template <class Range2, class ContractOperation>
+	MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperation<T,Operation,MAOps...> >
+	contract(const ContractOperation& cop, const std::string& indexName) const;
+
+	
+	template <class Range2, class ContractOperation>
+	MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperation<T,Operation,MAOps...> >
+	contract(const ContractOperation& cop, const std::string& indexName,
+		 const typename Range2::IndexType& begin,
+		 const typename Range2::IndexType& end) const;
+	*/
 	
 	template <class MAOp2>
 	auto operator+(const MAOp2& sec) -> decltype(operator()(std::plus<T>(), sec));
@@ -279,21 +296,23 @@ namespace MultiArrayTools
 	OBT mArgs; // include first arg also here !!!
     };
 
-    template <typename T, class ContractOperation, class Range, class MAOp>
-    class MultiArrayContraction : public MultiArrayOperation<T,ContractOperation,MAOp>
+    template <typename T, class ContractOperation, class Range, class... MAOps>
+    class MultiArrayContraction : public MultiArrayOperation<T,ContractOperation,MAOps...>
     {
     public:
 	typedef MultiArrayOperationBase<T> MAOB;
-	typedef MultiArrayOperation<T,ContractOperation,MAOp> MAO;
+	typedef MultiArrayOperation<T,ContractOperation,MAOps...> MAO;
 
-	MultiArrayContraction(const ContractOperation& cop, const MAOp& mao,
-			      const typename Range::IndexType& runIndex);
+	MultiArrayContraction(const ContractOperation& cop,
+			      const typename Range::IndexType& runIndex,
+			      const MAOps&... mao);
 
 	
-	MultiArrayContraction(const ContractOperation& cop, const MAOp& mao,
+	MultiArrayContraction(const ContractOperation& cop,
 			      const typename Range::IndexType& runIndex,
 			      const typename Range::IndexType& beginIndex,
-			      const typename Range::IndexType& endIndex);
+			      const typename Range::IndexType& endIndex,
+			      const MAOps&... mao);
 	
 	virtual const T& get() const override;
 
diff --git a/src/unit_test.cc b/src/unit_test.cc
index f1c7298..ae64fec 100644
--- a/src/unit_test.cc
+++ b/src/unit_test.cc
@@ -26,7 +26,7 @@ namespace {
 	    return res += a;
 	}
 	
-	void endOp(T& res) const {}
+	void endOp() const {}
 
 	void reset() const
 	{