Col_functions.cc
No OneTemporary
Actions

Size

176 KB

Subscribers

None

Col_functions.cc
View Options

	/*
	* Col_functions.cc
	* Contains definitions of functions and operators used for treating the color structure
	* Created on: Jul 8, 2010
	* Author: malin
	*/

	#include "Col_functions.h"


	namespace ColorFull {

	//using namespace std;
	using std::tr1::shared_ptr;

	/////////////////////// DEFINING OPERATORS /////////////////////

	// Defining + operator to be able to erase elements at place
	std::list<int>::iterator operator+( std::list<int>::iterator x, int n ) {

	while(n>0){
	x++;
	n--;
	}
	while(n<0){
	x--;
	n++;
	}
	return x;
	}

	std::list<int>::iterator operator-( std::list<int>::iterator x, int n ) {

	while(n>0){
	x--;
	n--;
	}
	while(n<0){
	x++;
	n++;
	}
	return x;
	}

	// Defining + operator to be able to erase elements at place
	std::list< Quark_line >::iterator operator+( std::list < Quark_line >::iterator x, int n ){

	while(n>0){
	x++;
	n--;
	}
	while(n<0){
	x--;
	n++;
	}
	return x;
	}

	col_str::iterator operator-( col_str::iterator x, int n ){

	while(n>0){
	x--;
	n--;
	}
	while(n<0){
	x++;
	n++;
	}
	return x;
	}


	// Define the operator << for vector of int
	std::ostream& operator<<( std::ostream& out, std::vector<int> vec ){
	int max=vec.size();
	if(max==0) out <<"{}";
	else{
	out <<"{";
	for (int i=0; i<max-1; i++){
	out << vec.at(i) << ",";
	}
	out << vec.at(max-1) <<"}";
	}
	return out;
	}

	// Define the operator << for cvec
	std::ostream& operator<<( std::ostream& out, const cvec & cv ) {

	out << "{";
	// Loop over entries
	for (uint i = 0; i < cv.size(); i++) {
	// Print element
	std::cout.width(6);
	std::ostringstream outstr;
	outstr << cv.at(i);
	out << outstr.str();
	// If not last element print ","
	if (i < (cv.size() -1 )) out << ", ";
	}
	out << "}";
	return out;
	}

	// Define the operator << for dvec
	std::ostream& operator<<(std::ostream& out, const dvec & dv) {

	out << "{";
	// Loop over entries
	for (uint i = 0; i < dv.size(); i++) {
	// Print element
	std::cout.width(6);
	std::ostringstream outstr;
	outstr << dv.at(i);
	out << outstr.str();
	// If not last element print ","
	if (i < (dv.size() -1 )) out << ", ";
	}
	out << "}";
	return out;
	}

	// Define the operator << for cmatr
	std::ostream& operator<<( std::ostream& out, const cmatr & cm ){

	out <<"{" << std::endl;
	// Loop over rows
	for(uint i=0; i< cm.size(); i++ ){
	out <<"{";
	// Loop over columns
	for(uint j=0; j< cm.at(i).size(); j++ ){
	// Print element
	std::cout.width( 6 );
	std::ostringstream outstr;
	outstr << cm.at(i).at(j);
	// If not last element print ","
	if (j<cm.at(i).size()-1 ) outstr << ",";
	out << outstr.str();
	//out << Poly_m.at(i).at(j) << "";
	}
	out <<"}";
	// If not last row, print ","
	if (i<cm.at(i).size()-1 ) out << ",";
	out << std::endl;
	}
	out <<"}" << std::endl;
	return out;
	}

	// Define the operator << for dmatr
	std::ostream& operator<<( std::ostream& out, const dmatr & matr ){

	//std::cout.precision(10);
	out <<"{" << std::endl;
	// Loop over rows
	for(uint i=0; i< matr.size(); i++ ){
	out <<"{";
	// Loop over columns
	for(uint j=0; j< matr.at(i).size(); j++ ){
	// Print element
	std::cout.width( 16 );
	std::ostringstream outstr;
	outstr.width( 20 );
	outstr.precision(16);
	// If the result is larger than accuracy print it out
	if( std::abs( matr.at(i).at(j) )> accuracy ){
	outstr << std::fixed << matr.at(i).at(j);
	}
	// otherwise is should probably be 0
	else
	outstr << std::fixed << 0;
	// If not last element print ","
	if (j<matr.at(i).size()-1 ) outstr << ",";
	out << outstr.str();
	//out << Poly_m.at(i).at(j) << "";
	}
	out <<"}";
	// If not last row, print ","
	if (i<matr.at(i).size()-1 ) out << ",";
	out << std::endl;
	}
	out <<"}" <<std::endl;
	return out;
	}


	// Define the operator << for dvec
	std::ostream& operator<<( std::ostream& out, std::pair<int, int> pair ) {

	out << "(";
	out << pair.first;
	out << ", ";
	out << pair.second;
	out << ")";

	return out;
	}

	/////////////////////// DEFINING FUNCKTIONS ////////////////////


	// Function for emitting a gluon from a Col_str
	Col_amp Col_functions::emit_gluon( const Col_str & in_Col_str, int emitter, int g_new ) const{

	// Locate the emitter in the Col_str
	std::vector <int> place=in_Col_str.find_parton(emitter);

	// Find what kind, q qbar or g, the emitter is
	std::string kind=in_Col_str.find_kind(emitter);

	// Defining Col_str to return (two needed in case of g)
	Col_str out_Col_str1=in_Col_str;
	Col_str out_Col_str2=in_Col_str;
	// Defining Col_amp to return
	Col_amp out_Col_amp;

	// If the emitter is a quark
	if( kind =="q" ){
	// Add new parton index at second first place in relevant quark_line
	out_Col_str1.insert( place.at(0), place.at(1)+1, g_new );
	out_Col_amp.ca.push_back( out_Col_str1 );
	}
	// If the emitter is an anti-quark
	else if( kind =="qbar" ){
	// Add new parton after the qbar
	out_Col_str1.insert( place.at(0), place.at(1), g_new );
	// Change sign
	// TODO check change of sign
	out_Col_str1.Poly=out_Col_str1.Poly*(-1);
	// Making a Col_amp to return
	out_Col_amp.ca.push_back( out_Col_str1 );
	}
	// If the emitter is a g
	else if( kind =="g" ){
	// If the emitter is a gluon a new gluon line should be inserted in two
	// different places, before and after the emitter

	// check sign, currently -sign when inserting the gluon before
	// Change sign Poly for first term
	Monomial Mon_tmp;
	Mon_tmp.int_part=-1;
	out_Col_str1.Poly=out_Col_str1.Poly*Mon_tmp;

	// Inserting the gluon before
	out_Col_str1.insert( place.at(0), place.at(1), g_new );
	// Inserting the gluon after
	out_Col_str2.insert( place.at(0), place.at(1)+1, g_new );
	// Appending result to out_Col_amp
	out_Col_amp.ca.push_back( out_Col_str1 );
	out_Col_amp.ca.push_back( out_Col_str2 );
	// Normal order
	out_Col_amp.normal_order();
	};
	out_Col_amp.simplify();

	return out_Col_amp;
	}

	// Function for emitting a gluon from a Col_amp
	Col_amp Col_functions::emit_gluon( const Col_amp & Ca_in, int emitter, int g_new ) const{

	Col_amp Ca_out;

	// Special case of gluons only, and only 2 gluons
	// As the anti-cyclic order is implicit, the third gluon should
	// only be added in one place.
	// For >2 gluons this will be automatic as the starting state contains
	// only one of the orders, for example (1,2,3) but not (1,3,2)
	// If there is only one Col_str, and that Cs has only one Ql,
	// and that Ql has only 2 partons and there are only gluons
	if(Ca_in.ca.size()==1 && Ca_in.ca.at(0).cs.size()==1 && Ca_in.ca.at(0).cs.at(0).ql.size()==2 && Ca_in.gluons_only()){

	// Ca_in.ca.size()==1 && Ca_in.gluons_only()Ca_in.ca.at(0).cs.size()==1 && Ca_in.ca.at(0).cs.at(0).ql.size()==2 &&
	// The new Col_amp has only one Col_str, which is typically (1,2,3)
	// This should give the old Col_str
	Ca_out=Ca_in;
	// Then the third gluon should be appended to the only ql in the only cs
	Ca_out.ca.at(0).cs.at(0).append(g_new);
	//cout << "emit_gluon(Ca..): Special case of 2 gluons, to return " << Ca_out.ca.at(0).cs.at(0);
	// check sign

	}
	else{ // Emit from each Col_str, and append to new Col_amp
	for( uint m=0; m< Ca_in.ca.size(); m++ ){
	// Emit from the Col_str
	Col_amp part_m=emit_gluon(Ca_in.ca.at(m), emitter, g_new);
	// Append result
	Ca_out.append(part_m.ca);
	}
	}
	//cout << "emit_gluon(Ca..): after emission" << endl << Ca_out << endl;


	return Ca_out;
	}
	/*
	// Function for splitting an open Quark_line into two Quark_lines
	// j is the locations of the gluons to be removed in the split
	Col_str Col_functions::split_after( Quark_line& Ql, int j ) const{

	if(!Ql.open){
	std::cerr <<"split_after: expects an open quark_line" << std::endl;
	std::cerr.flush();
	throw -1;
	}
	if( j==0 or j== static_cast<int>(Ql.ql.size()-1) ){
	std::cerr <<"split_after: the argument must be a gluon place, was " << j << std::endl;
	std::cerr.flush();
	throw -1;
	}

	Quark_line Ql1=Ql;
	Quark_line Ql2=Ql;

	// In first part erase everything after j
	quark_line::iterator it1=Ql1.ql.begin()+j+1;
	quark_line::iterator it2=Ql1.ql.end();
	Ql1.ql.erase(it1, it2);

	// In second part, erase everyting before j
	it1=Ql2.ql.begin();
	it2=Ql2.ql.begin()+j+1;
	Ql2.ql.erase( it1, it2 );


	Col_str Cs;

	// Construct Col_str to return
	Cs.cs.push_back(Ql1);
	Cs.cs.push_back(Ql2);

	return Cs;
	}
	*/
	/*
	// Returns a Quark_line with the the elements before j in a Quark_line
	Quark_line Col_functions::before( Quark_line Ql, int j ) const{
	//cout << "before( Quark_line Ql, int j ): the size was " << static_cast<int>(Ql.ql.size());
	//cout << "for the Quark_line " << Ql << endl;
	if( j<0 or j> static_cast<int>(Ql.ql.size()-1) ){
	std::cerr << "before( Quark_line Ql, int j ): the size was " << static_cast<int>(Ql.ql.size());
	std::cerr << "for the Quark_line " << Ql << std::endl;
	std::cerr.flush();
	std::cerr <<"before: the argument must be >=0 and inside vector, was " << j << std::endl;
	std::cerr.flush();
	throw -1;
	}

	// In first part erase everything after and includingj
	quark_line::iterator it1=Ql.ql.begin()+j;
	quark_line::iterator it2=Ql.ql.end();
	Ql.ql.erase(it1, it2);

	return Ql;
	}

	// Returns the elements after j in a Quark_line
	Quark_line Col_functions::after( Quark_line & Ql, int j ) const{

	if(j<0 or j> static_cast<int>(Ql.ql.size()-1)){
	std::cerr <<"after: the argument must be >=0 and inside vector, was " << j << std::endl;
	std::cerr.flush();
	throw -1;
	}

	// In first part erase everything after and including j
	quark_line::iterator it1=Ql.ql.begin();
	quark_line::iterator it2=Ql.ql.begin()+j+1;
	Ql.ql.erase(it1, it2);

	return Ql;
	}
	*/

	// Function for exchanging a gluon between the partons p1 and p2
	// See my general color structure paper
	Col_amp Col_functions::exchange_gluon( const Col_str & Cs, int p1, int p2 ) const{

	Col_str Cs_copy=Cs;

	// Find out kind and location of partons
	std::string kind1 = Cs_copy.find_kind(p1);
	std::string kind2 = Cs_copy.find_kind(p2);
	std::vector<int> place1 = Cs_copy.find_parton(p1);
	std::vector<int> place2 = Cs_copy.find_parton(p2);



	// To contain result
	Col_amp Ca;

	// Make sure the col_f in multiplying participating Ql's are 0
	Polynomial Poly1; // For comparing Poly is 1, all powers are 0
	if (!(Poly1 == Cs_copy.cs.at(place1.at(0)).Poly)) {
	// Move factor of Quark_line to Col_str
	Cs_copy.Poly = Cs_copy.Poly * Cs_copy.cs.at(place1.at(0)).Poly;
	Cs_copy.cs.at(place1.at(0)).Poly = Poly1;
	}
	if (!(Cs_copy.cs.at(place2.at(0)).Poly == Poly1)) {
	// Move factor of Quark_line to Col_str
	Cs_copy.Poly = Cs_copy.Poly * Cs_copy.cs.at(place2.at(0)).Poly;
	Cs_copy.cs.at(place2.at(0)).Poly = Poly1;
	}

	// If the qs are the same, the result is just CF*old for q and qbar, or (1/2)
	// sign
	if ( place1 == place2 ) {
	Monomial mon;
	// For qq or qbar qbar
	if (kind1 == "q" or kind1 == "qbar")
	mon.pow_CF = 1;
	// for gluons
	else
	mon.pow_Nc = 1;
	// There is a factor TR from the contracted gluon,
	// and we have 2 terms which survives in the large Nc limit
	// the factor -1 can be seen by noting that once the gluon
	// is inserted to the left and once to the right in the surviving terms
	// TODO check again!
	mon.pow_TR+=1;
	mon.int_part*=-2;

	Cs_copy.Poly = Cs_copy.Poly * mon;
	Ca.ca.push_back( Cs_copy );
	return Ca;
	}

	// Exchange between qq or qbar qbar
	if ((kind1 == "q" && kind2 == "q") or (kind1 == "qbar" && kind2 == "qbar")) {

	Col_str Cs2 = Cs_copy;

	// The q's normally sit on different quark_lines

	// The leading N-tems has the q's exchanged
	Cs_copy.cs.at(place1.at(0)).ql.at(place1.at((1))) = p2;
	Cs_copy.cs.at(place2.at(0)).ql.at(place2.at((1))) = p1;
	// Multiply with TR 1/2
	Monomial Mon;
	//Mon.pow_2 = -1; changed 121218
	Mon.pow_TR = 1;

	Cs_copy.Poly = Cs_copy.Poly * Mon;

	// The suppressed term is a copy of the old term,
	// but multiplied with -TR /Nc
	//Mon.pow_2 = -1; changed 121218
	Mon.pow_TR = 1;
	Mon.int_part = -1;
	Mon.pow_Nc = -1;
	Cs2.Poly = Cs2.Poly * Mon;

	Ca.ca.push_back(Cs_copy);
	Ca.ca.push_back(Cs2);
	}

	// Exchange between q qbar
	if ((kind1 == "q" && kind2 == "qbar") or (kind1 == "qbar" && kind2 == "q")) {

	std::vector<int> q_place;
	std::vector<int> qbar_place;
	int the_q;
	int the_qbar;

	if ((kind1 == "q" && kind2 == "qbar")) {
	q_place = place1;
	qbar_place = place2;
	the_q = p1;
	the_qbar = p2;
	} else if ((kind1 == "qbar" && kind2 == "q")) {
	q_place = place2;
	qbar_place = place1;
	the_q = p2;
	the_qbar = p1;
	}

	Col_str Cs2 = Cs_copy;

	// If q and qbar are not part of same ql
	if (q_place.at(0) != qbar_place.at(0)) {
	// The non-suppressed term, obtained by erasing q qbar and joining rest
	Cs_copy.erase(q_place);
	Cs_copy.erase(qbar_place);

	// In the ql of the qbar, append the ql of the q, then erase the ql of the q
	//cout << "Before append " << Cs_copy << endl;
	Cs_copy.cs.at(qbar_place.at(0)).append(Cs_copy.cs.at(q_place.at(0)).ql);
	//cout << "After append " << Cs_copy << endl;
	// and erase the ql of the q
	Cs_copy.erase(q_place.at(0));

	}
	// If q and qbar sit on same ql
	else {
	//std::cout << "q qbar part of same ql " << std::endl;
	// The non-suppressed term, obtained by erasing q qbar and joining rest to a ring
	Cs_copy.erase(qbar_place);
	Cs_copy.erase(q_place);
	Cs_copy.cs.at(q_place.at(0)).open = false;
	// If only one quark remains the term should be 0
	//if ( Cs_copy.cs.at( q_place.at(0) ).ql.size()==1 ) Cs_copy.clear();
	}
	// Construct Quark_line={q,qbar}, and append it
	Quark_line Ql;
	Ql.open = true;
	Ql.ql.push_back(the_q);
	Ql.ql.push_back(the_qbar);
	Cs_copy.cs.push_back(Ql);

	// Multiply with 1/2
	Monomial Mon;
	// Mon.pow_2 = -1; // Changed 121218
	Mon.pow_TR = 1;

	Cs_copy.Poly = Cs_copy.Poly * Mon;
	//std::cout << "Non-suppressed term " << Cs_copy;

	// The suppressed term is a copy of the old term,
	// but multiplied with -TR / Nc
	//Mon.pow_2 = -1;// Changed 121218
	Mon.pow_TR=1;
	Mon.int_part = -1;
	Mon.pow_Nc = -1;
	Cs2.Poly = Cs2.Poly * Mon;

	// Construct the resulting color_amplitude and return
	// Make sure there is not only one or 0 gluons in a closed ql
	// Append Cs to ca if >=2 gluons
	// if open, append
	if (Cs_copy.cs.at(q_place.at(0)).open)
	Ca.ca.push_back(Cs_copy);
	// if closed but >= 2 gluons, append
	else if (Cs_copy.cs.at(q_place.at(0)).ql.size() >= 2 && !Cs_copy.cs.at(
	q_place.at(0)).open)
	Ca.ca.push_back(Cs_copy);
	// If 1 gluon, the term is 0 and is not appended

	// If 0 gluons the empty colsed ql is equivalent to multiplying with Nc, obtained after simplification
	else if (Cs_copy.cs.at(q_place.at(0)).ql.empty()
	&& !Cs_copy.cs.at(q_place.at(0)).open) {
	// return 2nd term * (1-Nc^2)
	//Monomial Mon_tmp; // old incorrect code
	//Mon_tmp.pow_Nc = 2;
	//Cs2.Poly = Cs2.Poly + Mon_tmp * Cs2.Poly;
	Ca.ca.push_back(Cs_copy);
	Ca.ca.push_back(Cs2);

	Ca.simplify();
	//cout << "Ca in exchange_gluon" << Ca;
	return Ca;
	}
	Ca.ca.push_back(Cs2);
	}

	// Exchange between q g
	if ((kind1 == "q" && kind2 == "g") or (kind1 == "g" && kind2 == "q")) {

	std::vector<int> q_place;
	std::vector<int> g_place;
	int the_q=0;
	int the_g=0;

	if (kind1 == "q" && kind2 == "g") {
	q_place = place1;
	g_place = place2;
	the_q = p1;
	the_g = p2;
	} else if (kind1 == "g" && kind2 == "q") {
	q_place = place2;
	g_place = place1;
	the_q = p2;
	the_g = p1;
	}

	// For containing Quark_lines
	Quark_line Ql11;
	Quark_line Ql12;
	Quark_line Ql21;
	Quark_line Ql22;

	Col_str Cs1 = Cs_copy;
	Col_str Cs2 = Cs_copy;

	// If the q and the g are not part of same ql
	if (g_place.at(0) != q_place.at(0)) {

	// Take split the Ql with the gluon after the gluon
	Quark_line Qlg_start = Cs_copy.cs.at(g_place.at(0)).before( g_place.at(1) );

	Quark_line Qlg_end = Cs_copy.cs.at(g_place.at(0)).after( g_place.at(1) );
	Quark_line Qlq_end = Cs_copy.cs.at(q_place.at(0)).after( q_place.at(1) );

	// First part, + sign
	// First part, first Ql
	Ql11 = Qlg_end;
	Ql11.prepend(the_g);
	Ql11.prepend(the_q);

	// First part, second Ql
	Ql12 = Qlg_start;
	Ql12.append(Qlq_end.ql);

	// If the ql of the gluon was closed account for this by gluing qls together
	if (!Qlg_start.open) {
	Ql11.append(Ql12.ql);
	Ql11.open = true;
	};

	// Second part, first Ql
	Ql21 = Qlg_end;
	Ql21.prepend(the_q);

	// Second part, second Ql
	Ql22 = Qlg_start;
	Ql22.append(the_g);
	Ql22.append(Qlq_end.ql);

	// If the ql of the gluon was closed account for this by gluing qls together
	if (!Qlg_start.open) {
	Ql21.append(Ql22.ql);
	Ql21.open = true;
	};

	// Replace Ql's in the Col_str
	Cs1.cs.at(q_place.at(0)) = Ql11;
	Cs2.cs.at(q_place.at(0)) = Ql21;

	// There is only one ql if the ql of the gluon was closed
	if (Qlg_start.open) {
	Cs1.cs.at(g_place.at(0)) = Ql12;
	Cs2.cs.at(g_place.at(0)) = Ql22;
	} else { // Erase 2nd ql
	Cs1.erase(g_place.at(0));
	Cs2.erase(g_place.at(0));
	}

	}
	// If q and g are part of same ql
	else {

	Quark_line Ql1 = Cs_copy.cs.at(g_place.at(0)).before( g_place.at(1) );
	Ql1.ql.erase(Ql1.ql.begin());
	Quark_line Ql2 = Cs_copy.cs.at(g_place.at(0)).after( g_place.at(1));

	// cout << "Ql1 " << Ql1 << endl;
	// cout << "Ql2 " << Ql2 << endl;
	// First term, First ql
	Ql11 = Ql1;
	Ql11.open = false;

	// First term, 2nd ql
	Ql12 = Ql2;
	Ql12.prepend(the_g);
	Ql12.prepend(the_q);

	// 2nd term, 1st ql
	Ql21 = Ql1;
	Ql21.append(the_g);
	Ql21.open = false;

	// 2nd term 2nd ql
	Ql22 = Ql2;
	Ql22.prepend(the_q);

	// Replace Ql's in the Col_str
	Cs1.cs.at(q_place.at(0)) = Ql11;
	Cs1.cs.insert(Cs1.cs.begin() + q_place.at(0) + 1, Ql12);

	// Replace Ql's in the Col_str
	Cs2.cs.at(q_place.at(0)) = Ql21;
	Cs2.cs.insert(Cs2.cs.begin() + q_place.at(0) + 1, Ql22);

	}

	// Multiply first part with 1/2
	Monomial Mon_tmp;
	//Mon_tmp.pow_2 = -1; // Changed 121218
	Mon_tmp.pow_TR = 1;

	Cs1.Poly = Cs1.Poly * Mon_tmp;

	// Multiply second part with -1/2
	Monomial Mon_tmp2;
	//Mon_tmp2.pow_2 = -1;// Changed 121218
	Mon_tmp2.pow_TR = 1;
	Mon_tmp2.int_part = -1;
	Cs2.Poly = Cs2.Poly * Mon_tmp2;


	// Add to result
	Ca.ca.push_back(Cs1);
	Ca.ca.push_back(Cs2);
	}// end echange between q g


	// Exchange between qbar g (q and g in my paper)
	if ((kind1 == "qbar" && kind2 == "g") or (kind1 == "g" && kind2 == "qbar")) {

	//cout <<"qbar g \n";
	std::vector<int> qbar_place;
	std::vector<int> g_place;
	int the_qbar=0;
	int the_g=0;

	if (kind1 == "qbar" && kind2 == "g") {
	qbar_place = place1;
	g_place = place2;
	the_qbar = p1;
	the_g = p2;
	} else if (kind1 == "g" && kind2 == "qbar") {
	qbar_place = place2;
	g_place = place1;
	the_qbar = p2;
	the_g = p1;
	}

	// For containing Quark_lines
	Quark_line Ql11, Ql12, Ql21, Ql22;

	// For containing resulting Col_str's
	Col_str Cs1 = Cs_copy;
	Col_str Cs2 = Cs_copy;

	// If ql and g on different quark_lines
	if (g_place.at(0) != qbar_place.at(0)) {
	// Split ql's into parts
	Quark_line Qlg_start = Cs_copy.cs.at(g_place.at(0)).before( g_place.at(1)) ;
	Quark_line Qlg_end = Cs_copy.cs.at(g_place.at(0)).after( g_place.at(1) );
	Quark_line Qlqbar_start = Cs_copy.cs.at( qbar_place.at(0) ).before(qbar_place.at(1));

	// First part, + sign, first Ql
	Ql11 = Qlqbar_start;
	Ql11.append(the_g);
	Ql11.append(Qlg_end.ql);

	// First part, + sign, second Ql
	Ql12 = Qlg_start;
	Ql12.append(the_qbar);

	// If the ql of the gluon was closed account for this by gluing qls together
	if (!Qlg_start.open) {
	Ql11.append(Ql12.ql);
	Ql11.open = true;
	//cout << "Ql11" << Ql11 << endl;
	std::cout.flush();
	};

	// Second part, - sign, first Ql
	Ql21 = Qlqbar_start;
	Ql21.append(Qlg_end.ql);

	// Second part, - sign, second Ql
	Ql22 = Qlg_start;
	Ql22.append(the_g);
	Ql22.append(the_qbar);

	// If the ql of the gluon was closed account for this by gluing qls together
	if (!Qlg_start.open) {
	Ql21.append(Ql22.ql);
	Ql21.open = true;
	//cout << "Ql21" << Ql21 << endl;
	std::cout.flush();
	};

	// Replace Ql's in the Col_str
	Cs1 = Cs_copy;
	Cs2 = Cs_copy;
	if (Qlg_start.open) {
	Cs1.cs.at(qbar_place.at(0)) = Ql11;
	Cs1.cs.at(g_place.at(0)) = Ql12;
	Cs2.cs.at(qbar_place.at(0)) = Ql21;
	Cs2.cs.at(g_place.at(0)) = Ql22;
	} else { // If the ql of the g was closed, there is only one ql, erase 2nd ql
	Cs1.cs.at(qbar_place.at(0)) = Ql11;
	//cout << "Cs1 " << Cs1 << endl;
	Cs1.erase(g_place.at(0));

	Cs2.cs.at(qbar_place.at(0)) = Ql21;
	Cs2.erase(g_place.at(0));
	//cout << "Cs1 " << Cs1 << endl;
	//cout << "Cs2 " << Cs2 << endl;
	}
	}
	// If qbar and g are part of same quark_line
	else {

	Quark_line Ql1 = Cs_copy.cs.at(g_place.at(0)).before(g_place.at(1) );
	Quark_line Ql2 = Cs_copy.cs.at(g_place.at(0)).after( g_place.at(1) );
	Ql2.ql.erase(Ql2.ql.end() - 1);

	// First term, First ql
	Ql11 = Ql1;
	Ql11.append(the_qbar);

	// First term, 2nd ql
	Ql12 = Ql2;
	Ql12.prepend(the_g);
	Ql12.open = false;

	// 2nd term, 1st ql
	Ql21 = Ql1;
	Ql21.append(the_g);
	Ql21.append(the_qbar);

	// 2nd term 2nd ql
	Ql22 = Ql2;
	Ql22.open = false;

	// Replace Ql's in the Col_str
	Cs1 = Cs_copy;
	Cs1.cs.at(qbar_place.at(0)) = Ql11;
	Cs1.cs.insert(Cs1.cs.begin() + qbar_place.at(0) + 1, Ql12);

	// Replace Ql's in the Col_str
	Cs2 = Cs_copy;
	Cs2.cs.at(qbar_place.at(0)) = Ql21;
	Cs2.cs.insert(Cs2.cs.begin() + qbar_place.at(0) + 1, Ql22);
	}
	// Multiply with 1/2
	Monomial Mon_tmp;
	//Mon_tmp.pow_2 = -1; changed 121218
	Mon_tmp.pow_TR = 1;

	Cs1.Poly = Cs1.Poly * Mon_tmp;

	// Multiply with -1/2
	Monomial Mon_tmp2;
	//Mon_tmp2.pow_2 = -1; changed 121218
	Mon_tmp2.pow_TR = 1;
	Mon_tmp2.int_part = -1;
	Cs2.Poly = Cs2.Poly * Mon_tmp2;

	Ca.ca.push_back(Cs1);
	Ca.ca.push_back(Cs2);

	}

	// Exchange between g g
	if (kind1 == "g" && kind2 == "g") {

	// For containing resulting Col_str's
	Col_str Cs1 = Cs_copy;
	Col_str Cs2 = Cs_copy;
	Col_str Cs3 = Cs_copy;
	Col_str Cs4 = Cs_copy;

	// For containing Quark_lines
	Quark_line Ql11, Ql12, Ql21, Ql22, Ql31, Ql32, Ql41, Ql42;

	// If g's on different Ql's
	if (place1.at(0) != place2.at(0)) {

	// Split both ql's into parts
	Quark_line Ql1_start = Cs_copy.cs.at(place1.at(0)).before( place1.at(1) );
	Quark_line Ql1_end = Cs_copy.cs.at(place1.at(0)).after( place1.at(1) );
	Quark_line Ql2_start = Cs_copy.cs.at(place2.at(0)).before( place2.at(1) );
	Quark_line Ql2_end = Cs_copy.cs.at(place2.at(0)).after( place2.at(1) );

	// First term, both quarks on 2nd ql, First ql
	Ql11 = Ql1_start;
	Ql11.append(Ql2_end.ql);
	// First term, 2nd ql
	Ql12 = Ql2_start;
	Ql12.append(p2);
	Ql12.append(p1);
	Ql12.append(Ql1_end.ql);

	// 2nd term, g1 on ql1, g2 on ql2, 1st ql
	Ql21 = Ql1_start;
	Ql21.append(p1);
	Ql21.append(Ql2_end.ql);
	// 2nd term, 2nd ql
	Ql22 = Ql2_start;
	Ql22.append(p2);
	Ql22.append(Ql1_end.ql);

	// 3rd term, g2 on ql1, g1 on ql2,1st ql
	Ql31 = Ql1_start;
	Ql31.append(p2);
	Ql31.append(Ql2_end.ql);
	// 3rd term, 2nd ql
	Ql32 = Ql2_start;
	Ql32.append(p1);
	Ql32.append(Ql1_end.ql);

	// last term, both quarks on 1st ql,1st ql
	Ql41 = Ql1_start;
	Ql41.append(p1);
	Ql41.append(p2);
	Ql41.append(Ql2_end.ql);

	// last term, 2nd ql
	Ql42 = Ql2_start;
	Ql42.append(Ql1_end.ql);

	// Replace Ql's in the Col_str
	// If both Ql's are open
	if (Cs_copy.cs.at(place1.at(0)).open && Cs_copy.cs.at(place2.at(0)).open) {
	Cs1.cs.at(place1.at(0)) = Ql11;
	Cs1.cs.at(place2.at(0)) = Ql12;

	Cs2.cs.at(place1.at(0)) = Ql21;
	Cs2.cs.at(place2.at(0)) = Ql22;

	Cs3.cs.at(place1.at(0)) = Ql31;
	Cs3.cs.at(place2.at(0)) = Ql32;

	Cs4.cs.at(place1.at(0)) = Ql41;
	Cs4.cs.at(place2.at(0)) = Ql42;
	}
	// If string where closed glue together ends
	// If first ql closed
	else if (!Cs_copy.cs.at(place1.at(0)).open) {
	Ql12.append(Ql11.ql);
	Ql22.append(Ql21.ql);
	Ql32.append(Ql31.ql);
	Ql42.append(Ql41.ql);

	//The result should be open/closed if the ql2 was open/closed
	Ql12.open = Cs_copy.cs.at(place2.at(0)).open;
	Ql22.open = Cs_copy.cs.at(place2.at(0)).open;
	Ql32.open = Cs_copy.cs.at(place2.at(0)).open;
	Ql42.open = Cs_copy.cs.at(place2.at(0)).open;

	//The reslut contaied int Qlx2 sould be used for the Color structure
	// Replace Ql's in the Col_str
	Cs1.cs.at(place2.at(0)) = Ql12;
	Cs2.cs.at(place2.at(0)) = Ql22;
	Cs3.cs.at(place2.at(0)) = Ql32;
	Cs4.cs.at(place2.at(0)) = Ql42;

	Cs1.erase(place1.at(0));
	Cs2.erase(place1.at(0));
	Cs3.erase(place1.at(0));
	Cs4.erase(place1.at(0));
	}
	// if only 2nd ql closed
	if (!Cs_copy.cs.at(place2.at(0)).open && Cs_copy.cs.at(place1.at(0)).open) {
	Ql11.append(Ql12.ql);
	Ql21.append(Ql22.ql);
	Ql31.append(Ql32.ql);
	Ql41.append(Ql42.ql);

	// The result should be open
	Ql11.open = true;
	Ql21.open = true;
	Ql21.open = true;
	Ql21.open = true;

	//The result contained int Qlx1 should be used for the Color structure
	// Replace Ql's in the Col_str
	Cs1.cs.at(place2.at(0)) = Ql11;
	Cs2.cs.at(place2.at(0)) = Ql21;
	Cs3.cs.at(place2.at(0)) = Ql31;
	Cs4.cs.at(place2.at(0)) = Ql41;

	Cs1.erase(place1.at(0));
	Cs2.erase(place1.at(0));
	Cs3.erase(place1.at(0));
	Cs4.erase(place1.at(0));

	}
	}
	// If g's on same ql's
	if (place1.at(0) == place2.at(0)) {

	//cout << "g's on same ql \n";
	// Split ql into parts

	Quark_line Ql_start = Cs_copy.cs.at(place1.at(0)).before( std::min(place2.at(
	1), place1.at(1)));
	// the part between the g's
	Quark_line Ql_between = Cs_copy.cs.at(place1.at(0)).before(std::max(
	place2.at(1), place1.at(1)));
	Ql_between = Ql_between.after( std::min(place2.at(1), place1.at(1)) );

	// the end after 2nd g
	Quark_line Ql_end = Cs_copy.cs.at(place1.at(0)).after( std::max(place2.at(1),
	place1.at(1)));

	//std::cout << "QL_start " << Ql_start << std::endl;
	//std::cout << "Ql_between " << Ql_between << std::endl;
	//std::cout << "Ql_end " << Ql_end << std::endl;

	// added 130710
	// we have to know what g is first to put them back in right order
	int first_g, last_g;
	if( place1.at(1) < place2.at(1) ) {
	first_g=p1;
	last_g=p2;
	}
	else {
	first_g=p2;
	last_g=p1;
	}


	// First term, both gluons inbetween
	Ql11 = Ql_start;
	Ql11.append(Ql_end.ql);
	// First term, part between quarks
	Ql12 = Ql_between;
	//Ql12.append(p2);
	//Ql12.append(p1); changed 130710
	Ql12.append( last_g );
	Ql12.append( first_g );
	Ql12.open = false;

	//std::cout << "Ql11 " << Ql11 << std::endl;
	//std::cout << "Ql12 " << Ql12 << std::endl;


	// 2nd term, g1 after first part, g2 after 2nd
	Ql21 = Ql_start;
	//Ql21.append(p1); changed 130710
	Ql21.append( first_g );
	Ql21.append(Ql_end.ql);
	// First term, part between quarks
	Ql22 = Ql_between;
	//Ql22.append(p2); changed 130710
	Ql22.append( last_g );
	Ql22.open = false;

	//std::cout << "Ql21 " << Ql21 << std::endl;
	//std::cout << "Ql22 " << Ql22 << std::endl;

	// 3rd term, g2 after first part, g1 after 2nd
	Ql31 = Ql_start;
	//Ql31.append(p2); changed 130710
	Ql31.append( last_g);
	Ql31.append(Ql_end.ql);
	// First term, part between quarks
	Ql32 = Ql_between;
	//Ql32.append(p1);changed 130710
	Ql32.append( first_g );
	Ql32.open = false;

	//std::cout << "Ql31 " << Ql31 << std::endl;
	//std::cout << "Ql32 " << Ql32 << std::endl;

	// 4th term, g2 after first part, g1 after 2nd
	Ql41 = Ql_start;
	Ql41.append( first_g );
	Ql41.append( last_g );
	Ql41.append(Ql_end.ql);
	// First term, part between quarks
	Ql42 = Ql_between;
	Ql42.open = false;

	//std::cout << "Ql41 " << Ql41 << std::endl;
	//std::cout << "Ql42 " << Ql42 << std::endl;


	Cs1.cs.at(place1.at(0)) = Ql11;
	Cs1.cs.insert(Cs1.cs.begin() + place1.at(0) + 1, Ql12);
	Cs2.cs.at(place1.at(0)) = Ql21;
	Cs2.cs.insert(Cs2.cs.begin() + place1.at(0) + 1, Ql22);
	Cs3.cs.at(place1.at(0)) = Ql31;
	Cs3.cs.insert(Cs3.cs.begin() + place1.at(0) + 1, Ql32);
	Cs4.cs.at(place1.at(0)) = Ql41;
	Cs4.cs.insert(Cs4.cs.begin() + place1.at(0) + 1, Ql42);

	}

	// Multiplying with 1/2 where appropriate
	Monomial Mon_tmp;
	//Mon_tmp.pow_2 = -1;
	Mon_tmp.pow_TR = 1;

	Cs2.Poly = Cs2.Poly * Mon_tmp;
	Cs3.Poly = Cs3.Poly * Mon_tmp;

	// Multiplying with -1/2 where appropriate
	//Mon_tmp.pow_2 = -1; changed 121218
	Mon_tmp.pow_TR = 1;

	Mon_tmp.int_part = -1;
	Cs1.Poly = Cs1.Poly * Mon_tmp;
	Cs4.Poly = Cs4.Poly * Mon_tmp;

	Ca.ca.push_back(Cs1);
	Ca.ca.push_back(Cs2);
	Ca.ca.push_back(Cs3);
	Ca.ca.push_back(Cs4);

	}

	// Simplify Col_amp
	Ca.simplify();

	return Ca;

	}


	// Function for exchanging a gluon between two partons p1 and p2
	Col_amp Col_functions::exchange_gluon( const Col_amp & Ca, int p1, int p2 ) const{

	// Final col_amp to return
	Col_amp Ca_out;

	// Exchange in each Col_str, and append to new Col_amp
	for(uint m=0; m< Ca.ca.size(); m++ ){
	// Exchange in Col_str m
	Col_amp part_m=exchange_gluon( Ca.ca.at(m), p1, p2 );
	// Append result
	Ca_out.append(part_m.ca);
	}

	return Ca_out;
	}

	// Exchange gluons between all possible parton pairs
	Col_amp Col_functions::exchange_gluon( Col_amp & Ca, int n_p ) const{

	Col_amp Ca_out;

	// Loop over partons
	for(int p1=1; p1<= n_p; p1++){
	for(int p2=1; p2<= n_p; p2++){

	Ca_out=Ca_out+ exchange_gluon(Ca, p1, p2);

	}
	}
	Ca_out.simplify();
	return Ca_out;
	}


	// Function simulating virtual gluon exchane, evolving the color structure
	// Evolution performed n_steps time with size alpha
	Col_amp Col_functions::step_evolv( Col_amp & Ca, int n_p, int n_steps, double alpha ) const{

	// Let A->A+alpha*ANew, repeat this N_step times
	for(int i=1; i<=n_steps; i++){
	// The amplitude resulting after exchange, part to add
	Col_amp Ca_new=exchange_gluon(Ca, n_p);

	// The evolved amplitude, A->A+alpha*ANew
	Ca=Ca+Ca_new*alpha;
	Ca.simplify();
	}
	return Ca;
	}

	// Function for contracting quarks
	// between two color structures Cs1 and Cs2
	Col_str Col_functions::contract_quarks( Col_str Cs1, Col_str Cs2 ) const{
	//cout << "contract_quarks(Cs1,Cs2): has the arguments " << Cs1 << Cs2 << endl;
	// A vector for containing place of next quark to remove
	// in Cs1 and in conjugate( Col_str Cs2)
	std::vector<int> q_place;
	std::vector<int> q_place2;

	// The conjugate of Cs2
	//Col_str conj_Cs2 = conjugate(Cs2);

	// The conjugate of Cs1
	Col_str conj_Cs1 = Cs1;
	conj_Cs1.conjugate();

	// The total color structure
	//Col_str Cs = Cs1*conj_Cs2;
	// Changed to conjugate first arg, 130524
	Col_str Cs = conj_Cs1*Cs2;

	//cout << "contract_quarks(Cs1, Cs2): Cs1 Cs2^* " << Cs << endl;

	// Count how many quarks should be contracted
	int n_q = Cs.n_quark();

	// As long as there are quark_lines left to contract
	while (n_q > 0) {
	// Find first quark in Cs1 by looping over Quark_lines
	for (int i = 0; (n_q>0 && i < static_cast<int>(Cs.cs.size()) ); i++) {
	//cout << "contract_quarks(Cs1, Cs2): Looking at place " << i << " in " <<Cs << endl;
	std::cout.flush();
	// Check if the quark-line is open, in which case it has a q
	if (Cs.cs.at(i).open) {
	//cout << "contract_quarks(Cs1, Cs2): Found open ql at " << i
	// << " " << Cs.cs.at(i) << endl;
	// The first quark is located and has position
	q_place.clear();
	q_place.push_back(i);
	q_place.push_back(0);
	// and number
	int q = Cs.at(q_place.at(0), q_place.at(1));
	//cout << "contract_quarks(Cs1, Cs2): Found q to remove " << q << " at place " << q_place << endl;

	// Locate same quark a second time
	// Loop over Quark_lines
	q_place2.clear();
	int i2 = i + 1; // Quark_line of second occurrence
	while (q_place2.empty()) { // As long as quark not found a second time
	if (Cs.cs.at(i2).at(Cs.cs.at(i2).ql.size() - 1) == q) {// If quark found, store place
	q_place2.push_back(i2);
	q_place2.push_back(Cs.cs.at(i2).ql.size() - 1);
	}
	i2++;
	}
	if (q_place2.empty()) {
	std::cerr << "contract_quarks(Cs1, Cs2): Found q " << q
	<< " only once in " << Cs << std::endl;

	}

	//cout << "contract_quarks(Cs1, Cs2): Also found at place "
	// << q_place2 << endl;

	// Prepare new Quark_line
	// to be inserted at the place of found open Quark_line
	Quark_line new_Quark_line;
	Quark_line part2_new_Quark_line;
	// The first part of the new Quark_line should be the Quark_line
	// containing q in the conjugate
	new_Quark_line = Cs.cs.at(q_place2.at(0));

	// Erasing q in the end
	new_Quark_line.ql.erase(--new_Quark_line.ql.end());
	//cout << "contract_quarks(Cs1, Cs2): First part after removing contracted quark, new_Quark_line, from conj_Cs2 " << new_Quark_line <<"\n";
	part2_new_Quark_line = Cs.cs.at(q_place.at(0));

	// Erasing the q in the beginning
	part2_new_Quark_line.ql.erase(part2_new_Quark_line.ql.begin());
	//cout << "contract_quarks(Cs1, Cs2): Second part " << part2_new_Quark_line <<"\n";

	new_Quark_line.append(part2_new_Quark_line.ql);
	//cout << "contract_quarks(Cs1, Cs2): The new quark_line is " << new_Quark_line << "\n";

	// If the first q index and the last q(bar) index in the new
	// Quark_line is the same (and the Quark_line is "open"), the indices
	// should be removed and the Quark_line should be closed
	if (new_Quark_line.ql.at(0) == new_Quark_line.ql.at(
	new_Quark_line.ql.size() - 1) && new_Quark_line.open) {
	//cout <<"First = Last \n";
	// The string is closed
	new_Quark_line.open = false;
	// Remove last and first index
	new_Quark_line.ql.erase(--new_Quark_line.ql.end(),
	new_Quark_line.ql.end());
	// cout <<"The new quark_line is " << new_Quark_line <<"\n";
	new_Quark_line.ql.erase(new_Quark_line.ql.begin());
	//cout <<"contract_quarks(Cs1, Cs2): The new quark_line is " << new_Quark_line <<"\n";
	}

	// Inserting new Quark_line in the place of the old
	Cs.cs.at(i) = new_Quark_line;
	//cout << "contract_quarks(Cs1, Cs2): The Cs after insertion is " << Cs << "\n";
	// Remove quark_line with q in Cs
	Cs.cs.erase((Cs.cs.begin() + q_place2.at(0)));
	//cout << "contract_quarks(Cs1, Cs2): after removal " << Cs << "\n";
	i=-1; // reset to keep looking from the beginning in the new Cs (i will be increased to 0)
	}// end of if (open)

	n_q = Cs.n_quark();
	//cout << "contract_quarks(Cs1, Cs2): The new number of quarks is " << n_q << endl;

	} // end of for, loop over quark_lines

	}
	//cout << "contract_quarks(Cs1, Cs2): To return " << Cs << endl;
	std::cout.flush();
	return Cs;

	}

	// Function for contracting quarks
	// between two Col_amps.
	Col_amp Col_functions::contract_quarks( Col_amp Ca1, Col_amp Ca2 ) const{
	Col_amp Ca_res;

	// Make sure the Col_strs are not empty "[]"=1, as all indices contracted
	Ca1.remove_empty_Col_strs();
	Ca2.remove_empty_Col_strs();

	// Loop over Col_strs, and contract quarks between all possible combinations
	// Loop over Col_str's in Ca1
	for(uint m1=0; m1 < Ca1.ca.size(); m1++ ){
	// Loop over Col_strs in Ca2
	for(uint m2=0; m2 < Ca2.ca.size(); m2++ ){
	Ca_res.ca.push_back( contract_quarks( Ca1.ca.at(m1), Ca2.ca.at(m2)) );
	}
	}

	/*
	// Col_amps should not have uncontracted quarks and a Polynomial at the same time
	if( !Ca1.Scalar.size() && !Ca2.Scalar.empty() ){
	cerr << "contract_quarks(Ca1, Ca1): Col_amp with uncontracted quarks had Polynomial " << Ca1.Scalar << " and " << Ca2.Scalar << endl;
	throw -1;

	}
	// Add Polynomials
	Ca_res.Scalar=Ca1.Scalar+Ca2.Scalar;
	*/
	//cout << "contract_quarks(Ca1, Ca2): To return " << Ca_res << endl;
	return Ca_res;
	}

	// Function to contract neighboring gluons in a Quark_line starting at j
	Quark_line Col_functions::contract_neighboring_gluons( Quark_line & Ql, int j ) const {
	//cout << "contract_neighboring_gluons( Ql, j ): incoming " << Ql <<" j " << j << endl;
	// If asked for -1st element, check last
	if(j==-1) j=Ql.ql.size()-1;
	// If asked for last+1 element, check first
	if( j== static_cast<int>(Ql.ql.size()-1) ) j=0;
	// Keep contracting gluons as long as:
	// there are at least two gluons to contract
	// neighbors are equal
	// and j is not the last parton in the ql
	while( j< static_cast<int>(Ql.ql.size()-1) && (Ql.ql.size()) >=2 && (Ql.ql.at(j))==Ql.ql.at(j+1) ){

	quark_line::iterator it1=Ql.ql.begin()+j;
	quark_line::iterator it2=Ql.ql.begin()+j+2;
	// Removing neighboring gluons
	Ql.ql.erase(it1,it2);

	// This should increase the power of Cf
	Monomial Mon_tmp;
	Mon_tmp.pow_CF=1;
	Ql.Poly=Ql.Poly*Mon_tmp;
	if(j>2) j=j-2;
	// cout << "contract_neighboring_gluons: removed one pair, new Ql" << Ql << "and new j is" << j << endl;
	}
	// If j is the last parton and the Quark_line is closed and first=last
	// and there are at least 2 gluons
	while( static_cast<int>(Ql.ql.size())>=2 && j== static_cast<int>(Ql.ql.size()-1) && !Ql.open && Ql.ql.at(0)==Ql.ql.at( Ql.ql.size()-1) ){
	// cout << "contract_neighboring_gluons( Ql, j ): entering second while " << endl;

	quark_line::iterator it1=Ql.ql.begin();
	quark_line::iterator it2=Ql.ql.end()-1;
	Ql.ql.erase(it2);
	Ql.ql.erase(it1);

	//cout << "contract_neighboring_gluons( Ql, j ): removed last and first, new Ql" << Ql << "and new j is" << j << endl;
	//cout.flush();

	// This should increase the power of Cf
	Monomial Mon_tmp;
	Mon_tmp.pow_CF=1;
	Ql.Poly=Ql.Poly*Mon_tmp;


	// If now, all gluons are removed, increase Nc an open the ql (if not already open)
	if( Ql.ql.empty() && !Ql.open ){
	Monomial Mon_tmp;
	Mon_tmp.pow_Nc=1;
	Ql.Poly=Ql.Poly*Mon_tmp;
	Ql.open=true;
	}

	j=j-2;
	}

	//cout << "contract_neighboring_gluons( Ql, j ): to return" << Ql << endl;
	return Ql;
	}


	// Function to contract neighboring gluons in a Quark_line starting at place 0
	Quark_line Col_functions::contract_neighboring_gluons( Quark_line & Ql ) const{

	//cout << "contract_neighboring_gluons( Ql ): Incoming " << Ql << endl;
	// Counting powers of CF,
	// Should be increases once for every gluon that is contracted
	// initially put to dummy value
	int CF_pow_old=-9999;
	int CF_pow_new=0;

	// Keep looping over elements as long as CF_pow changes
	while( CF_pow_old!= CF_pow_new ){
	CF_pow_old=CF_pow_new;

	// Loop over elements in the Quark_line, starting at first parton
	// and ending with last
	for (uint j=0; j < Ql.ql.size(); j++){
	Ql=contract_neighboring_gluons(Ql, j);
	}
	// To see if CF_pow changed, arbitrarily compare 0th element in Polynomial (if existing)
	// (all pow_CF changes)
	if( !Ql.Poly.empty() )
	CF_pow_new=Ql.Poly.at(0).pow_CF; // put CF_pow_new to current pow
	}
	//cout << "contract_neighboring_gluons( Ql ): To return " << Ql << endl;
	return Ql;
	}


	// Function to contract next to neighboring gluons in a Quark_line
	// starting at place j (so gluon j and j+2)
	// Also looks for new neighbors
	Quark_line Col_functions::contract_next_neighboring_gluons( Quark_line & Ql, int j) const{
	//cout << "contract_next_neighboring_gluons (Ql, j): incoming " << Ql << " int, "<< j<< "\n";

	// If asked for -1st element, check last
	if(j==-1) j=Ql.ql.size()-1;
	// If asked for -2nd element, check second last
	if(j==-2) j=Ql.ql.size()-2;

	// Keep track of changes in number of gluons
	// add dummy number '1' to run at least once
	int old_size=Ql.ql.size()+1;

	// Loop over index in quark_line as long as the size keep changing,
	// or at least once
	while( old_size!= static_cast<int>(Ql.ql.size()) && (Ql.ql.size())>=4){
	old_size= (Ql.ql.size());
	// There are three cases:
	// the normal case, and, if the Quark_line is closed,
	// 2nd last =first (0th) and last =2nd (1st)
	bool normal_case=(Ql.ql.size()>=4 && static_cast<uint>(j)< Ql.ql.size()-2 && Ql.ql.at(j)==Ql.ql.at(j+2));
	bool second_last_case=( Ql.ql.size()>=4 && !Ql.open && static_cast<uint>(j)==Ql.ql.size()-2 && Ql.ql.at(j)==Ql.ql.at(0) );
	bool last_case=(Ql.ql.size()>=4 && !Ql.open && static_cast<uint>(j)==Ql.ql.size()-1 && Ql.ql.at(j)==Ql.ql.at(1));
	// See if next to neighbors are equal
	// The point of having a while is to check again if one pair was found
	while(normal_case or second_last_case or last_case) {
	//cout << "Found next to neighbor " << j << " " << Ql.ql.at(j) <<"\n";

	// Erase the gluons (at places depending on size)
	quark_line::iterator it1;
	quark_line::iterator it2;
	if(normal_case){
	it1=Ql.ql.begin()+j;
	it2=Ql.ql.begin()+j+2;
	//cout << "Normal case \n";
	}
	else if(second_last_case){
	//recently changed
	it1=Ql.ql.begin();
	it2=Ql.ql.end()-2;
	//cout << "Second last case \n";
	}
	else if(last_case){
	it1=Ql.ql.begin()+1;
	it2=Ql.ql.end()-1;
	// Erase right most gluon first
	//cout << "Last case \n";

	}
	// Erase right most gluon first
	Ql.ql.erase(it2);
	Ql.ql.erase(it1);
	// The surviving term comes with -TR/(Nc)
	Monomial Mon_tmp;
	// Mon_tmp.pow_2=-1;
	Mon_tmp.pow_TR = 1;

	Mon_tmp.int_part=-1;
	Mon_tmp.pow_Nc=-1;
	//cout <<"Ql " << Ql <<endl;
	Ql.Poly=Ql.Poly*Mon_tmp;
	//cout << Ql << endl;

	// Check if new neighbors are equal
	uint old_size2=Ql.ql.size();
	// Check if inbetween gluon is now neighbor with itself to the right
	Ql=contract_neighboring_gluons(Ql, j);
	// Check if inbetween gluon is now neighbor with itself to the left
	Ql=contract_neighboring_gluons(Ql, j-1);
	// if a neighboring pair was removed
	if( (Ql.ql.size()) != old_size2 ) j=j-2;

	// Two gluons have been erased so to keep looking forward in the
	// Quark_line j must not be increased
	j--;
	// Why did it not work to change type of Ql.ql.size() using static_cast<int>??
	normal_case=( Ql.ql.size()>=4 && static_cast<uint>(j)< Ql.ql.size()-2 && Ql.ql.at(j)==Ql.ql.at(j+2));
	if(normal_case){
	//cout << "j "<< j << ", Ql.ql.size() "<< Ql.ql.size() << endl;
	}
	second_last_case=( Ql.ql.size()>=4 && !Ql.open && j== static_cast<int>(Ql.ql.size())-2 && Ql.ql.at(j)==Ql.ql.at(0));
	last_case=( Ql.ql.size()>=4 && !Ql.open && j== static_cast<int>(Ql.ql.size())-1 && Ql.ql.at(j)==Ql.ql.at(1));
	}
	}

	//cout << "contract_next_neighboring_gluons (Ql, j): to return " << Ql << "\n";
	return Ql;
	}


	// Function to contract next to neighboring gluons in a Quark_line
	Quark_line Col_functions::contract_next_neighboring_gluons( Quark_line & Ql) const{
	//cout << "contract_next_neighboring_gluons( Ql ): Incoming " << Ql << "\n";

	// Start with looking for neighbors
	Ql=contract_neighboring_gluons(Ql);
	// Keep track of changes in number of gluons
	// add dummy number '1' to run at least once
	int old_size=Ql.ql.size()+1;

	// Loop over index in quark-line as long as the size keep changing,
	// or at least once
	while( old_size!= static_cast<int>(Ql.ql.size()) && (Ql.ql.size())>=4){
	// For seeing if size has changed
	old_size=Ql.ql.size();
	for (int j=0; j< static_cast<int>(Ql.ql.size()); j++){
	// Contract next to neighbors starting at place j
	Ql=contract_next_neighboring_gluons(Ql, j);
	}

	}
	//cout << "contract_next_neighboring_gluons( Ql ): To return " << Ql << "\n";
	return Ql;
	}

	// Contract up to next to neighboring gluons in each Quark_line in the Col_str
	Col_str Col_functions::contract_next_neighboring_gluons( Col_str & Cs ) const{
	//cout << "contract_next_neighboring_gluons( Cs ): Incoming " << Cs <<endl;

	// Loop over Quark_lines and remove neighboring and next to neighboring
	// gluon indices
	for( uint i=0; i < Cs.cs.size(); i++ ){

	// Create a Quark_line to use as argument
	Quark_line Ql= Cs.at(i);

	Ql=contract_next_neighboring_gluons( Ql );

	// Adding powers to get total power of 2, Nc and CF
	// Collect these powers in Cs
	Cs.Poly=Cs.Poly*Ql.Poly;
	// and put powers in Ql to 0
	Ql.Poly.clear();

	Cs.cs.insert(Cs.cs.begin() +i, Ql);

	// Remove old version of Quark_line (now at place i+1)
	Cs.cs.erase( Cs.cs.begin() +i + 1 );
	}
	// Remove 1 and 0-rings, simplify Poly and normal order
	Cs.simplify();

	//cout << "contract_next_neighboring_gluons( Cs ): To return " << Cs <<endl;
	return Cs;
	}


	// Contract up to next to neighboring gluons in each Quark_line in the Col_str in each Col_amp
	Col_amp Col_functions::contract_next_neighboring_gluons( Col_amp & Ca ) const{
	// Loop over Col_str's
	for (uint m = 0; m < Ca.ca.size(); m++) {
	Ca.ca.at(m) = contract_next_neighboring_gluons(Ca.ca.at(m));
	}

	return Ca;
	}

	// Function for contracting gluon indices in a two-ring in a Col_str
	Col_str Col_functions::contract_2_rings( Col_str & Cs ) const{
	//cout << "contract_2_rings(Col_str): Called with argument " << Cs << endl;
	// Contract gluons from 2-ring, this gives only one term
	// For storing place of two-ring
	std::vector<int> place1;
	// For storing place of gluon with same index as second gluon in two-ring
	// i.e. the gluon to be replaced with the first index in the 2-ring
	std::vector<int> place2;
	// The second index in two-ring, to be removed
	int the_g;

	// Loop over Quark_lines to find a two-ring
	for (uint i = 0; i < Cs.cs.size(); i++) {
	place1.clear();
	place2.clear();
	// Search for 2-rings
	// A gluon 2-ring has length 2 and is closed
	if (Cs.cs.at(i).ql.size() == 2 && !Cs.cs.at(i).open) {

	// Save place of two-ring
	place1.push_back(i);
	// Pick second gluon (for no good reason)
	// this index should be replaced by first index
	// in the other place where it occurs, and the 2-ring should be removed
	place1.push_back(1);
	//cout << "Found 2-ring at "<< i << "\n";
	the_g = Cs.at(place1.at(0), place1.at(1));
	//cout << "The gluon to be reomved is " << the_g << " aat place " << place1;
	//std::cout.flush();
	// Now, in all Quark_lines, look for same the_g until found
	for (uint i2 = 0; i2 < Cs.cs.size(); i2++) {
	for (uint j2 = 0; j2 < Cs.cs.at(i2).ql.size(); j2++) {
	// If the_g was found at a DIFFERENT place
	if ((i2 != i or j2 != 1) && Cs.at(i2, j2) == the_g) {
	place2.push_back(i2);
	place2.push_back(j2);
	//cout << "Found same g, " << the_g << ", at place " << place2;
	}
	}
	}

	// Check that the gluon index was found again
	if (place2.empty()) {
	std::cerr << "contract_2_rings: Did not find index " << the_g
	<< " other than at place " << place1 << "\n";
	std::cerr.flush();
	throw -1;
	}
	// If the_g was found twice in the same ql
	if (place1.at(0) == place2.at(0)) {
	//cout << "Found gluon " << the_g << "twice in the quark_line "<< place1.at(0) << endl;

	// Keep the Monomial
	// Multiply with Nc Mon from the contraction
	Monomial Mon_tmp;
	Mon_tmp.pow_Nc = 1;
	Mon_tmp.pow_CF = 1;
	// Multiply the Poly of the Col_str with the Poly of the ql
	// and this color factor
	Cs.Poly = Cs.Poly * Cs.cs.at(place1.at(0)).Poly * Mon_tmp;

	// Erase the ql
	Cs.erase(place1.at(0));

	return Cs;
	}

	// That index should be changed to the index of the first gluon
	// in the 2-ring
	//cout << "The gluon to be removed is: " << Cs.at2( place2.at(0), place2.at(1) ) << "\n";
	Cs.cs.at(place2.at(0)).ql.at(place2.at(1))
	= Cs.at(place1.at(0), 0);
	//cout << "Replacement done with result: " << Cs << "\n";
	// The two ring should be removed
	Cs.cs.erase(Cs.cs.begin() + i);
	Monomial Mon_tmp;
	//Mon_tmp.pow_2 = -1;
	Mon_tmp.pow_TR = 1;

	Cs.Poly = Cs.Poly * Mon_tmp;
	i++;
	//cout << "2-ring removed with result: " << Cs << "\n";
	}
	}
	// One-rings may have been created
	Cs.remove_1_rings();

	//cout << "contract_2_rings(Col_str): To return " << Cs << endl;
	return Cs;
	}

	// Contract two-rings in each Quark_line in the Col_str in each Col_amp
	Col_amp Col_functions::contract_2_rings( Col_amp & Ca ) const{
	// Loop over Col_str's
	for (uint m = 0; m < Ca.ca.size(); m++) {
	Ca.ca.at(m) = contract_2_rings(Ca.ca.at(m));
	}
	return Ca;
	}

	// Function for splitting a closed Quark_line into two Quark_lines
	// j1 && j2 are the locations of the gluons to be removed in the split
	// May create 1-rings and 0-rings
	std::string Col_functions::spm_file_name( const Col_amp & Basis, bool leading ) const{

	// First construct filename
	std::ostringstream ss;
	std::string filename;
	ss << "ColorResults";
	ss << '/';
	ss << "SPM_q";
	int nq=Basis.at(0).n_quark();
	ss << nq;
	ss << "_g";
	int ng=Basis.at(0).n_gluon();
	ss << ng;
	if ( leading ) ss << "_l";
	if ( leading_CF ) ss << "_CFl";
	ss << "_Nc";
	ss << Nc;
	ss << ".dat";
	filename=ss.str();
	return filename;
	}

	// Function for contracting gluon indices within the same Quark_line
	// Checks only for ONE pair
	Col_amp Col_functions::contract_Quark_line_gluons( Quark_line & Ql ) const{
	//cout <<"contract_Ql_gluons(Ql): Incoming " << Ql << endl;

	// Check that all quark indices have been removed
	if (Ql.open){
	std::cerr
	<< "Col_functions::contract_Quark_line_gluons: all quark indices where not contracted in " << Ql << std::endl;
	}
	// There will be two terms to keep track of, store in Col_amp
	Col_amp Ca;

	//If the Quark_line was empty, return a Col_amp with info on Col_str
	if (Ql.ql.empty()) {
	Col_str Cs;
	// Move Polynomial to Col_str
	Cs.Poly = Ql.Poly;
	Ql.Poly.clear();
	Cs.cs.push_back(Ql);
	Ca.ca.push_back(Cs);
	return Ca;
	}

	// Too keep track if partner was found or not
	bool found_pair = false;

	// Take first term and look for partner
	for (uint j1 = 0; j1 < Ql.ql.size() - 1; j1++) {
	int the_g = Ql.ql.at(j1);
	// Look for same g
	for (uint j2 = j1 + 1; j2 < Ql.ql.size(); j2++) {
	// if same g found
	if (Ql.ql.at(j2) == the_g) {

	// Special case that the gluons are neighbors
	if (j2 == j1 + 1 or (!Ql.open && j1 == 0 && j2 == Ql.ql.size()
	- 1)) {
	Ql = contract_neighboring_gluons(Ql, j1);
	// Make a Col_str out of the Quark_line
	Col_str Cs;
	Cs.cs.push_back(Ql);
	Ca.ca.push_back(Cs);
	Ca.simplify();
	//cout << "contract_Ql_gluons(Ql): to return " << Ca << endl;
	return Ca;
	}

	// Special case that the gluons are next to neighbors
	if (j2 == j1 + 2 //normal case
	or (!Ql.open && j1 == 0 && j2 == Ql.ql.size()- 2) //first and second last
	or (!Ql.open && j1 == 1 && j2 == Ql.ql.size()- 1) //second and last
	) {
	//cout << "contract_Ql_gluons: Found next to neighbors in Ql"
	// << Ql << endl;
	if (j2 == j1 + 2)
	Ql = contract_next_neighboring_gluons(Ql, j1);
	if ((!Ql.open && j1 == 0 && j2 == Ql.ql.size() - 2))
	Ql = contract_next_neighboring_gluons(Ql, Ql.ql.size() - 2);
	if( (!Ql.open && j1 == 1 && j2 == Ql.ql.size()- 1) )
	Ql = contract_next_neighboring_gluons(Ql, Ql.ql.size() - 1);

	// Make a Col_str out of the Quark_line
	Col_str Cs;
	Cs.cs.push_back(Ql);
	Ca.ca.push_back(Cs);
	Ca.simplify();
	//cout << "contract_Ql_gluons(Ql): to return " << Ca << endl;
	return Ca;
	}

	// Make a Col_str out of the Quark_line
	Col_str Cs;

	// Move color factor to multiply Cs instead of Ql
	Cs.Poly = Ql.Poly;
	Ql.Poly.clear();
	Cs.cs.push_back(Ql);

	// Make two copies to store both terms
	Ca.ca.push_back(Cs);
	Ca.ca.push_back(Cs);

	// The first non-suppressed term, obtained by splitting the Quark_line
	//Ca.ca.at(0) = split_Quark_line(Ca.ca.at(0).cs.at(0), j1, j2);
	//Ca.ca.at(0) = Ca.ca.at(0).cs.at(0).split_Quark_line( j1, j2 );

	std::pair <Quark_line, Quark_line> Ql_parts=Ca.ca.at(0).cs.at(0).split_Quark_line( j1, j2 );
	Col_str Cs_Ql_parts;
	Cs_Ql_parts.cs.push_back(Ql_parts.first);
	Cs_Ql_parts.cs.push_back(Ql_parts.second);

	Ca.ca.at(0)=Cs_Ql_parts;

	// This is multiplying a factor TR
	Monomial Mon_tmp;
	Mon_tmp.pow_TR = 1;

	Ca.ca.at(0).Poly = Ca.ca.at(0).Poly * Mon_tmp;
	// The split can generate new (next to) neighbors
	//cout << "Looking for nn in" << Ca.at(0) << endl ;
	Ca.ca.at(0) = contract_next_neighboring_gluons(Ca.ca.at(0));

	//cout << "contract_Ql_gluons(Ql): First term + non changed" << Ca << endl ;

	// The second Nc suppressed term, obtained by just removing gluon
	// remove g
	quark_line::iterator it1 = Ca.ca.at(1).cs.at(0).ql.begin() + j1;
	quark_line::iterator it2 = Ca.ca.at(1).cs.at(0).ql.begin() + j2;
	Ca.ca.at(1).cs.at(0).ql.erase(it2);
	Ca.ca.at(1).cs.at(0).ql.erase(it1);
	// Multiply with -TR/(Nc)
	//Mon_tmp.pow_2 = -1;
	Mon_tmp.pow_TR = 1;
	Mon_tmp.pow_Nc = -1;
	Mon_tmp.int_part = -1;
	// cout <<" Ca.at(1).Poly old" << Ca.at(1).Poly <<endl;
	//cout <<" Mon_tmp " << Mon_tmp <<endl;
	Ca.ca.at(1).Poly = Ca.ca.at(1).Poly * Mon_tmp;

	//cout << "contract_Ql_gluons(Ql): both terms" << Ca << endl ;
	//cout <<" Ca.at(1).Poly new" << Ca.at(1).Poly <<endl;
	// Change sign
	// A pair was found, use for stoping
	found_pair = true;
	}
	if (found_pair)
	break; // Stop j2 loop
	}
	if (found_pair)
	break; // Stop j1 loop
	}

	// If a pair was never found return the info of the original Ql
	// but with the Polynomial moved to the Col_str
	if(!found_pair){
	Col_str Cs;
	Cs.Poly=Ql.Poly;
	Ql.Poly.clear();
	Cs.cs.push_back(Ql);
	Ca.ca.push_back(Cs);
	}

	// Remove possible 0 and 1-rings
	Ca.remove_1_rings();
	Ca.remove_0_rings();

	// cout <<"contract_Ql_gluons(Ql): To return ";
	//std::cout.flush();
	//cout << Ca << endl;
	//std::cout.flush();

	return Ca;
	}


	// Function for contracting gluon indices within the same Quark_line
	Col_amp Col_functions::contract_Quark_line_gluons( Col_amp & Ca ) const{
	//cout <<"contract_Ql_gluons(Ca): Incoming " << Ca << endl;
	//std::cout.flush();

	// If Col_amp empty, do nothing
	if (Ca.ca.empty())
	return Ca;

	// Loop over individual Col_str's
	for (uint m = 0; m < Ca.ca.size(); m++) {

	// In each Col_str, loop over individual Quark_lines
	for (uint i = 0; (m < Ca.ca.size() && i < Ca.ca.at(m).cs.size()); i++) {
	// the Ql to be contracted
	Quark_line Ql = Ca.ca.at(m).cs.at(i);
	//cout << "contract_Ql_gluons(Ca): Ql to be contracted for i= " << i << " " << Ql << endl;

	// Contract gluons in the Quark_line
	// If the ql had a Polynomial, this is contained in Ca_from_Ql
	Col_amp Ca_from_Ql;
	Ca_from_Ql = contract_Quark_line_gluons(Ql);
	//cout <<"contract_Ql_gluons(Ca): after contract_Ql_gluons( Ql ) in Quark_line " << i << " in Col_str m "<< m << " " << Ca_from_Ql << endl;

	// The whole Col_str has to be replaced with a linear
	// combination of different Col_str's where the Ql has changed to
	// the different options

	// The rest of the Col_str, without changed Ql
	// If the whole Col_str (at m) had a Polynomial, this is contained in
	// Cs_rest
	Col_str Cs_rest;
	Ca.ca.at(m).cs.erase( Ca.ca.at(m).cs.begin() + i);
	Cs_rest = Ca.ca.at(m);
	//cout <<"contract_Ql_gluons(Ca): Cs_rest " << Cs_rest << endl;

	// The rest of the Ca, without affected Cs
	Ca.erase(m);

	// Append Quark_line(s) to the rest of the Col_str,
	// Loop over Col_str's in Ca_from_Ql
	for (uint mq = 0; mq < Ca_from_Ql.ca.size(); mq++) {
	Col_str Cs_new = Cs_rest;
	// Append the col_str
	Cs_new.append(Ca_from_Ql.ca.at(mq).cs);
	//cout <<"contract_Ql_gluons(Ca): Cs_new from mq=" << mq << ", " << Cs_new << endl;
	// and multiply with the Polynomial
	Cs_new.Poly = Cs_new.Poly * Ca_from_Ql.ca.at(mq).Poly;

	//Append the new Col_str to the Ca
	Ca.ca.push_back(Cs_new);
	}
	// If the contraction gave a pure polynomial, this should also
	// give rise to a Col_str
	for (int mq = 0; mq < Ca_from_Ql.Scalar.size(); mq++) {
	// Unless the term in the Polynomial is 0
	if( Ca_from_Ql.Scalar.at(mq).int_part != 0 ){
	Col_str Cs_new = Cs_rest;
	// The Col_str contains no col_str,
	// but the old polynomial is multiplied with the Polynomial in Scalar
	Cs_new.Poly = Cs_new.Poly *Ca_from_Ql.Scalar.at(mq); // Problem when renaming Poly to Scalar?
	//cout <<"contract_Ql_gluons(Ca): Cs_new from mq=" << mq << ", " << Cs_new << endl;
	// Append the new Col_str to the Ca
	Ca.ca.push_back(Cs_new);
	}
	}
	//cout <<"contract_Ql_gluons(Ca): Ca after adding back terms " << Ca << endl;
	// If the contraction gave more than 1 term,
	// increase m too keep looking forward
	// TODO compensate for terms from Ca_from_Ql.Poly, m should change to keep
	// going forward
	if (Ca_from_Ql.ca.size() >= 2)
	m = m + Ca_from_Ql.ca.size() - 1;
	}
	}

	Ca.simplify();
	//cout << "contract_Ql_gluons(Ca): To return " << Ca << endl;
	return Ca;

	}

	// Function for contracting gluon indices within the same Quark_line in Col_str
	Col_amp Col_functions::contract_Quark_line_gluons( Col_str & Cs ) const{
	Col_amp Ca;
	Ca.ca.push_back(Cs);
	Ca= contract_Quark_line_gluons( Ca );
	return Ca;
	}


	// Contracts one gluon, the first g in first Quark_line
	// Treat only the case when the same g is found in different Ql
	Col_amp Col_functions::contract_a_gluon( Col_str & Cs ) const{
	//cout << "contract_a_gluon(Cs): Incoming " << Cs << endl;

	//If the Col_str is empty, or the first ql is empty return it as Col_amp
	if ( Cs.cs.empty() or Cs.cs.at(0).ql.empty() ) {
	Col_amp Ca;
	Cs.cs.clear();
	Ca.ca.push_back(Cs);
	return Ca;
	}

	// Pick first gluon
	int the_g = Cs.at(0, 0);
	std::vector<int> place;
	place.push_back(0);
	place.push_back(0);

	// For storing place of second gluon
	std::vector<int> place2;

	// Locate the same gluon
	for (uint i2 = 0; i2 < Cs.cs.size(); i2++) {
	for (uint j2 = 0; j2 < Cs.cs.at(i2).ql.size(); j2++) {
	//cout <<"i2 " << i2 <<" j2 " << j2 <<endl;
	// Make sure it is a different gluon
	if ((i2 != 0 or j2 != 0) && Cs.at(i2, j2) == the_g) {
	place2.push_back(i2);
	place2.push_back(j2);
	}
	}
	}

	// If the gluon was not found, something went wrong
	if (place2.empty()) {
	std::cerr << "contract_a_gluon: The gluon " << the_g
	<< " was only found once " << Cs;
	std::cerr.flush();
	throw -1;
	}

	// If the gluon was found in same Ql, use contract_Ql_gluons
	if (place.at(0) == place2.at(0)) {
	Col_amp Ca_out;
	Ca_out = contract_Quark_line_gluons(Cs);
	//cout << "contract_a_gluon(Cs): Gluon found twice in same Ql, returning contract_Ql_gluons( Cs ) " << Ca_out << endl;
	std::cout.flush();
	//cout << "contract_a_gluon(Cs): To return" << Ca_out;
	return Ca_out;
	}

	// If the removed gluon was part of a two-ring
	if( Cs.cs.at(place.at(0)).ql.size()==2 or Cs.cs.at(place2.at(0)).ql.size()==2){
	//cout << "contract_a_gluon(Cs): Found gluon in " << Cs << " is part of a two-ring "<< endl;
	Col_amp Ca_out;
	Cs=contract_2_rings(Cs);
	Ca_out.ca.push_back(Cs);
	//cout << "contract_a_gluon(Cs): After contracting two-rings " << Ca_out << endl;
	return Ca_out;

	}
	// The result is a sum of two terms, to contain these
	Col_str Cs1 = Cs;
	Col_str Cs2;

	// The Nc suppressed term, obtained by just removing the g
	// location of first and second g to remove
	quark_line::iterator it1 = Cs1.cs.at(place.at(0)).ql.begin();
	quark_line::iterator it2 = Cs1.cs.at(place2.at(0)).ql.begin() + place2.at(1);
	// Remove the gluon in both places
	Cs1.cs.at(place2.at(0)).ql.erase(it2);
	Cs1.cs.at(place.at(0)).ql.erase(it1);


	// The non-suppressed term, obtained by joining the ql's (where the g is removed)
	Cs2 = Cs1;
	// Construct the contracted ql by
	// 1, taking the first part of the 2nd ql, by erasing everything after the erased gluon
	Quark_line ql_first_part= Cs.cs.at(place2.at(0)).before( place2.at(1) );
	// 2, the inbetween part=first ql, the contracted gluon is already erased
	Quark_line middle_part=Cs2.cs.at(place.at(0));
	// 3, Last part equal to second part of 2nd ql
	Quark_line last_part= Cs.cs.at(place2.at(0)).after( place2.at(1) );
	// The new combined Ql, to replace the 2nd involved Qualk_line
	Quark_line new_Ql=ql_first_part;
	new_Ql.append( middle_part.ql );
	new_Ql.append( last_part.ql );
	// Copy polynomial info
	new_Ql.Poly=Cs.cs.at(place2.at(0)).Poly;
	// Replace the first involved Ql with the constructed Ql
	Cs2.cs.at(0)=new_Ql;
	//cout << "The contracted Ql is " << new_Ql<< endl;
	/*
	// Insert second term after first
	for (uint i = 0; i < Cs1.cs.at(place2.at(0)).ql.size(); i++) {
	Cs2.cs.at(0).ql.push_back(Cs2.cs.at(place2.at(0)).ql.at(i));
	}
	*/
	// erase moved ql
	col_str::iterator it = Cs2.cs.begin() + place2.at(0);
	Cs2.cs.erase(it);

	//cout <<"Cs2 " << Cs2 << endl;

	// Multiply with TR for non-suppressed term
	Monomial Mon_tmp;
	//Mon_tmp.pow_2 = -1;
	Mon_tmp.pow_TR = 1;
	Cs2.Poly = Cs2.Poly * Mon_tmp;

	// Multiply with -TR/Nc for suppressed term
	Mon_tmp.pow_Nc = -1;
	Mon_tmp.int_part = -1;
	Cs1.Poly = Cs1.Poly * Mon_tmp;

	// Look for simple simplifications in Cs1
	// This has to be done after Cs1 is used to construct Cs2
	// After the removal there may be new next to neighbors in the affected ql's
	// First affected Ql
	// look forward at nn
	Cs1.cs.at(place.at(0)) = contract_neighboring_gluons(Cs1.cs.at(place.at(0)),
	place.at(0));
	// look backward at nn
	Cs1.cs.at(place.at(0)) = contract_neighboring_gluons(Cs1.cs.at(place.at(0)),
	place.at(0) - 1);
	// look forward at next nei
	Cs1.cs.at(place.at(0)) = contract_next_neighboring_gluons(Cs1.cs.at(place.at(0)),
	place.at(0));
	// look backward at next nei
	Cs1.cs.at(place.at(0)) = contract_next_neighboring_gluons(Cs1.cs.at(place.at(0)),
	place.at(0) - 2);
	// Second affected ql
	// look forward at nn
	Cs1.cs.at(place2.at(0)) = contract_neighboring_gluons(Cs1.cs.at(place2.at(0)),
	place2.at(0));
	// look backward at nn
	Cs1.cs.at(place2.at(0)) = contract_neighboring_gluons(Cs1.cs.at(place2.at(0)),
	place2.at(0) - 1);
	// look forward at next nei
	Cs1.cs.at(place2.at(0)) = contract_next_neighboring_gluons(Cs1.cs.at(place2.at(0)),
	place2.at(0));
	// look backward at next nei
	Cs1.cs.at(place2.at(0)) = contract_next_neighboring_gluons(Cs1.cs.at(place2.at(0)),
	place2.at(0) - 2);

	Col_amp Ca_out;
	Ca_out.ca.push_back(Cs1);
	Ca_out.ca.push_back(Cs2);


	//cout << "contract_a_gluon(Cs): The resulting Ca before removing 1 and 0-rings " << Ca_out << endl;

	// remove rings with just one gluon (=0)
	Ca_out.remove_1_rings();
	// remove rings with no partons (=Nc, if closed)
	Ca_out.remove_0_rings();

	//cout << "contract_a_gluon(Cs): The resulting Ca before normal_ordering " << Ca_out << endl;
	std::cout.flush();
	//Ca_out.normal_order();
	//cout << "contract_a_gluon(Cs): The Ca to return is " << Ca_out << endl;
	std::cout.flush();
	return Ca_out;

	}

	// Contracts one gluon, the first g in first Quark_line (in each Col_str)
	Col_amp Col_functions::contract_all_gluons( Col_str & Cs ) const{
	//cout << "contract_all_gluons( Cs ): Incoming, " << Cs << endl;

	Col_amp Ca;
	Ca.ca.push_back(Cs);
	return contract_all_gluons(Ca);
	}

	Col_amp Col_functions::contract_a_gluon( Col_amp & Ca ) const{
	//cout << "on(Ca): Incoming " << Ca << endl;

	// If the Ca has only an empty Col_str or a Col_str with an empty Quark_line, do nothing
	if(Ca.ca.size()==1 && ( Ca.ca.at(0).cs.empty() or Ca.ca.at(0).cs.at(0).ql.empty())) return Ca;

	// Else collect contracted Col_amps in a variable Ca_res
	Col_amp Ca_res;
	Ca_res.Scalar=Ca.Scalar; // Move Polynomial

	// Make one contraction in each Col_str, move Col_str's when done to Ca_res
	while ( Ca.ca.size() > 0) {

	// Contract one gluon in Col_str
	Col_amp Ca_part=contract_a_gluon( Ca.ca.at(0) );
	//cout << "contract_a_gluon(Ca): Ca_part: " <<Ca <<endl;

	// ... and add resulting Col_amp to Ca_res
	Ca_res=Ca_res+Ca_part;
	// Erase old term (still at place 0)
	Ca.erase(0);
	}
	//cout << "contract_a_gluon(Ca): To return: " <<Ca <<endl;;
	return Ca_res;
	}

	// Function for contracting all gluon indices in a Col_str,
	// Assumes quarks already contracted
	Col_amp Col_functions::contract_all_gluons( Col_amp & Ca ) const{
	//cout << "contract_all_gluons( Ca ): Incoming, " << Ca << endl;

	// Make sure the Col_strs are not empty "[]"=1, as all indices contracted
	Ca.remove_empty_Col_strs();
	Ca.remove_empty_Col_strs();

	// Check that all quark indices have been removed
	if ( !Ca.gluons_only() ) {
	std::cerr << "contract_all_gluons(Cs): Error, all quark indices where not contracted in " << Ca << std::endl;
	std::cerr.flush();
	throw -1;
	}

	// First normal order and look for simple simplification
	Ca.simplify();

	Col_amp Ca_old; // For comparison
	// Contract remaining gluons as long as result keep changing
	// Start with contractions only giving one term
	while (Ca_old != Ca) {

	Ca_old = Ca;

	// Contract gluons from 2-ring as this gives only one term
	// and may result in new neighbors or next to neighbors
	Ca = contract_2_rings(Ca);
	//cout << "contract_all_gluons(Ca): contract_2_rings(Ca) done, with result Ca=" << Ca << endl;

	// Remove neighboring and next to neighboring gluon indices
	// in each quark_line as this gives only one term
	Ca = contract_next_neighboring_gluons(Ca);
	// cout << "contract_all_gluons(Ca): contract_next_neighboring_gluons done, with result Ca="
	// << Ca << endl;

	// Temporarely added to see if problem with empty Ql's disappeare
	//Ca.simplify();

	// Contract gluons within same quark_line
	// This will give >1 term, but at least rings will be shorter and shorter
	Ca = contract_Quark_line_gluons( Ca );

	// Recently introduces to maximize number of CF's
	Ca = contract_next_neighboring_gluons(Ca);


	// Make "arbitrary" gluon contraction when no simple or Ql contraction remains
	// In each Col_str, contract the first gluon in the first Ql
	//cout << endl << "contract_all_gluons(Ca): Ca before contract_a_gluon(Ca) " << Ca <<endl;
	Ca = contract_a_gluon(Ca);
	//cout.precision(20);
	//cout.setf(ios_base::fixed );
	//cout << endl << "contract_all_gluons(Ca): Ca after contract_a_gluon(Ca) " << Ca <<endl;
	//cout << "contract_all_gluons: Ca_old " << Ca_old << endl;

	// Look if same Col_str is represented more than once in the col_amp and simplify
	Ca.simplify();

	} //end of while (Ca keeps changing)

	//cout << "contract_all_gluons(Ca): Returning: " << Ca << endl;
	return Ca;
	}

	// Function for finding the leading power of Nc in a Polynomial
	int Col_functions::leading_Nc_pow( const Polynomial & Poly ) const{
	if( Poly.poly.empty()) {
	if( double_num(Poly)!=0 ) return 0;
	else return std::numeric_limits<int>::min();
	}
	int leading_pow=Poly.at(0).pow_Nc + Poly.at(0).pow_CF;

	// Loop over non-zero Monomials
	for (uint i=0; i<Poly.poly.size(); i++){
	if( Poly.at(i).pow_Nc + Poly.at(i).pow_CF > leading_pow && Poly.at(i).int_part!=0 )
	leading_pow=Poly.at(i).pow_Nc + Poly.at(i).pow_CF;
	}
	return leading_pow;
	}

	// Function for finding the leading power of Nc in a Polynomial
	int Col_functions::leading_Nc_pow( const Poly_vec & Pv ) const{
	int leading_pow=leading_Nc_pow( Pv.at(0) );

	// Loop over Polynomials
	for (uint p=0; p<Pv.size(); p++){
	if ( leading_Nc_pow( Pv.at(p) ) > leading_pow )
	leading_pow=leading_Nc_pow( Pv.at(p) );
	}

	return leading_pow;
	}

	// Function for finding the leading power of Nc among the Polynomials a vector
	// of Polynomial pointers points to
	int Col_functions::leading_Nc_pow( std::vector< std::tr1::shared_ptr<Polynomial> > Pvp) const{
	int leading_pow=leading_Nc_pow( *(Pvp.at(0)) );

	// Loop over Polynomials
	for (uint p=1; p<Pvp.size(); p++){
	if ( leading_Nc_pow( *(Pvp.at(p)) ) > leading_pow )
	leading_pow=leading_Nc_pow( *(Pvp.at(p)) );
	}

	return leading_pow;
	}

	// Takes the leading Nc terms of a Polynonmial, i.e. Monomials with highest power of Nc+CF
	// Replaces CF or not depending on the value of the global parameter leading_CF.
	Polynomial Col_functions::leading( const Polynomial & Poly, bool lcf ) const{

	// Candidate for highest power of Nc+CF
	int cand_pow = Poly.at(0).pow_Nc + Poly.at(0).pow_CF;

	// For containing the leading Nc Polynomial
	Polynomial leading_pol;
	leading_pol.push_back(Poly.at(0));

	// Looks for Monomials with highest Nc power
	// If more than one term, look for higher powers
	if( Poly.size()>1 ){// found bug for Polynomials with one term 13 07 12
	for (int k = 1; k < Poly.size(); k++) {

	// If power higher
	if ((Poly.at(k).pow_Nc + Poly.at(k).pow_CF) > cand_pow) {
	cand_pow = Poly.at(k).pow_Nc + Poly.at(k).pow_CF;
	leading_pol.clear();
	leading_pol = leading_pol * 0;
	leading_pol = leading_pol + Poly.at(k);
	}
	// If power equal
	else if (Poly.at(k).pow_Nc + Poly.at(k).pow_CF == cand_pow) {
	leading_pol = leading_pol + Poly.at(k);
	}
	}
	}
	// now leading_pol contains terms with maximal power of CF plus Nc

	// If the variable leading_CF is true, of if the function is called
	// with the argument lcf=true, CF should be replaced with TR*Nc.
	// If numerical evaluation is done after this, CF will be evaluated to TR Nc.
	// If CF is not replaced in this way, then, by numerical evaluation CF will
	// be replaced by CF(Nc), including the sub-leading term
	if ( leading_CF \|\| lcf ) {
	// Loop over terms and replace CF with the Nc -> infinity limit of CF =TR*Nc
	for (int i = 0; i < leading_pol.size(); i++) {
	int pow_CF = leading_pol.at(i).pow_CF;
	leading_pol.at(i).pow_Nc += pow_CF;
	leading_pol.at(i).pow_TR += pow_CF;
	leading_pol.at(i).pow_CF = 0;
	}
	}
	leading_pol.simplify();
	return leading_pol;
	}

	// Take the leading part of a Poly_vec
	Poly_vec Col_functions::leading( const Poly_vec & Pv ) const{

	// To contain the result
	Poly_vec Pv_res;

	// Loop over entries in vector (Polynomials), and take the leading part of each Polynomial
	// after this all Monomials in the SAME Polynomial have the same Nc power
	for (uint i = 0; i < Pv.size(); i++) {
	// Take the leading terms in each component
	Pv_res.push_back( leading( Pv.at(i) ) );
	}

	// Find the leading power
	int leading_pow = leading_Nc_pow(Pv_res);

	// Loop over entries and keep only those with maximal power
	for (uint p = 0; p < Pv_res.size(); p++) {
	if (leading_Nc_pow(Pv_res.at(p)) != leading_pow) {
	Pv_res.at(p) = 0 * Pv_res.at(p);
	Pv_res.at(p).simplify();
	}
	}

	return Pv_res;
	}

	// Take the leading part of a Poly_matr, the highest Nc powers in the whole matrix.
	Poly_matr Col_functions::leading( const Poly_matr & Pm ) const{

	// To contain the result
	Poly_matr Pm_res;

	// Loop over Poly_vecs, and take the leading part of each Poly_vec
	for (uint i = 0; i < Pm.size(); i++) {
	// Take the leading terms in each Poly_vec
	Poly_vec Pvl= leading(Pm.at(i));
	Pm_res.push_back( Pvl );
	}
	int cand_pow = leading_Nc_pow( Pm.at(0).pv );

	// Loop over Poly_vecs to locate highest power
	for (uint k = 0; k < Pm_res.size(); k++) {
	if (leading_Nc_pow( Pm_res.at(k).pv ) > cand_pow)
	cand_pow = leading_Nc_pow( Pm_res.at(k).pv );
	}

	// In all Polynomials, in all Poly_vec's,
	// put all elements which doesn't have highest power to 0
	for (uint k = 0; k < Pm_res.size(); k++) {
	// In each Poly_vec, loop over all Polynomials
	for (uint p = 0; p < Pm_res.at(k).size(); p++) {
	if (leading_Nc_pow(Pm_res.at(k).at(p)) != cand_pow) {
	//Pm_res.at(k).pv.at(p) = Pm_res.at(k).pv.at(p) * 0;
	Pm_res.at( k,p )=Pm_res.at(k).at(p) * 0;
	Pm_res.at(k,p).simplify();
	}
	}
	}

	return Pm_res;
	}


	// Take the leading part of a Poly_vec when given a pointer to it
	Poly_vec Col_functions::leading( std::vector<std::tr1::shared_ptr<Polynomial> > Pvp ) const{
	//cout << "leading(std::vector<std::vector<std::tr1::shared_ptr<Polynomial> > >Pm): Entering" << endl;
	Poly_vec Pv_res;

	// Loop over entries in vector (Poly_vecs), and take the leading part of each Polynomial
	// after this each Monomial has the same Nc+CF-power
	for (uint i = 0; i < Pvp.size(); i++) {
	// Take the leading terms in each component
	std::tr1::shared_ptr<Polynomial> the_pointer=Pvp.at(i);
	Pv_res.push_back( leading( (*the_pointer)) );
	}

	// Find the leading power
	int leading_pow = leading_Nc_pow(Pv_res);

	// Loop over entries and keep only those with maximal power
	for (uint p = 0; p < Pv_res.size(); p++) {
	if (leading_Nc_pow(Pv_res.at(p)) != leading_pow) {
	Pv_res.at(p) = 0 * Pv_res.at(p);
	Pv_res.at(p).simplify();
	}
	}

	return Pv_res;
	}

	// Take the leading part of a Poly_matr, the highest Nc powers in the whole matrix
	dmatr Col_functions::leading( std::vector< std::vector< std::tr1::shared_ptr<Polynomial> > > Ppm ) const {
	std::cout << "leading(std::vector<std::vector<std::tr1::shared_ptr<Polynomial> > >Pm): Finding leading terms in matrix..." << std::endl;

	// To contain the result as a matrix of double
	dmatr res;
	res.reserve( Ppm.size() );

	int cand_pow = leading_Nc_pow( Ppm.at(0) );

	// Loop over Poly_vecs to locate highest power
	for (uint k = 0; k < Ppm.size(); k++) {
	int the_pow=leading_Nc_pow(Ppm.at(k));
	if ( the_pow > cand_pow )
	cand_pow = the_pow;
	}
	//cout << "Found highest power " << cand_pow << endl;

	// In all Polynomials, in all Poly_vec's,
	// put all elements which doesn't have highest power to 0
	for (uint k = 0; k < Ppm.size(); k++) {
	//cout << "k is now " << k << endl;
	dvec dummy;
	res.push_back(dummy); // Not to run out of range
	// In each Poly_vec, loop over all Polynomials
	for (uint p = 0; p < Ppm.at(k).size(); p++) {
	Polynomial the_poly=*Ppm.at(k).at(p);
	leading_Nc_pow( the_poly );
	if ( leading_Nc_pow( the_poly ) != cand_pow ) {res.at(k).push_back(0);}
	else res.at(k).push_back( double_num(leading(the_poly)) );
	}
	}
	std::cout << " done." << std::endl;
	return res;
	}

	/// To keep only leading terms in a map.
	std::map< std::string, Polynomial > Col_functions::leading( std::map< std::string, Polynomial > mem_map ) const{

	//cout << "mem_map.size() " << mem_map.size() << endl;

	// To contain (string, leading(Polynomial))
	std::map< std::string, Polynomial > res;

	// Iterator for looping
	//std::map< string, tr1::shared_ptr<Polynomial> >::iterator iter=mem_map.begin();
	std::map< std::string, Polynomial >::iterator iter=mem_map.begin();


	// Find the highest power of Nc+CF
	//cout << "The poly is " << mem_map[ static_cast<string>(iter->first) ] << endl;
	//cout << "The poly is " << (iter->second) << endl;
	int pow_cand=leading_Nc_pow( (iter->second) );
	//cout << "pow_cand " << pow_cand;
	for( iter = mem_map.begin(); iter != mem_map.end(); ++iter ) {
	//cout << iter-> first <<endl;
	//cout.flush();
	Polynomial Poly= (iter->second);
	if( leading_Nc_pow(Poly) > pow_cand ) pow_cand=leading_Nc_pow(Poly);
	}
	//cout << "Winning pow_cand " << pow_cand;

	for( iter = mem_map.begin(); iter != mem_map.end(); ++iter ) {
	//cout << iter-> first <<endl;

	// Insert pair of string and the leading versions of the Polynomial
	Polynomial lead_Poly= leading((iter->second));
	if ( leading_Nc_pow(lead_Poly) != pow_cand ) lead_Poly=lead_Poly*0;

	//tr1::shared_ptr<Polynomial> ptr=shared_ptr<Polynomial> (new Polynomial(lead_Poly));
	res.insert( make_pair( iter->first, lead_Poly ) );
	} // After this only the leading part of each Polynomial contributes



	return res;
	}

	// Numerically evaluates a Monomial to a complex cnum number.
	cnum Col_functions::cnum_num( const Monomial & Mon ) const {

	cnum res=pow(Nc, Mon.pow_Nc)pow(CF, Mon.pow_CF)pow(TR, Mon.pow_TR)Mon.int_part Mon.cnum_part;

	return res;
	}

	// Numerically evaluates a Monomial to a double.
	double Col_functions::double_num( const Monomial & Mon ) const{

	double ratio =imag( cnum_num(Mon) )/real( cnum_num(Mon) );
	// Warn if the complex number has significant imaginary parts
	if( ratio > accuracy )
	std::cerr << "Warning keeping only real part of complex number, the ratio im/re was " << ratio << std::endl;

	return real( cnum_num(Mon) );
	}

	// Numerically evaluates a Polynomial for Nc=3
	cnum Col_functions::cnum_num( const Polynomial & Poly ) const {

	// An empty Polynomial has numerical value 0
	cnum res = 0;

	// Add contributions from Monomials
	for (int k = 0; k < Poly.size(); k++) {
	cnum part_k;
	part_k = pow(Nc, Poly.at(k).pow_Nc)
	*
	// The CF below is relevant only if CF has not been replaced in leading(Polynomial)
	pow(CF, Poly.at(k).pow_CF)
	*
	pow(TR, Poly.at(k).pow_TR) * Poly.at(k).int_part
	* Poly.at(k).cnum_part;
	res = res + part_k;
	}

	return res;
	}

	// Numerically evaluates a Polynomial.
	double Col_functions::double_num( const Polynomial & Poly ) const{

	double ratio =imag( cnum_num(Poly) )/real( cnum_num(Poly) );
	// Warn if the complex number has significant imaginary parts
	if( ratio > accuracy )
	std::cerr << "Warning keeping only real part of complex number, the ratio im/re was " << ratio << std::endl;

	return real( cnum_num(Poly) );
	}
	/*
	// Keeping only leading terms, numerically evaluate a Polynomial
	cnum Col_functions::leading_cnum_num(const Polynomial & Poly) const{
	cnum res=0;

	// Keep only leading terms in the polynomial
	// Poly=leading( Poly );
	Polynomial leading_Poly=leading( Poly );

	// Add contributions from remaining leading Monomials, putting
	for ( int k=0; k < leading_Poly.size(); k++ ){
	cnum part_k;
	part_k=pow( Nc, leading_Poly.at(k).pow_Nc )*
	pow( CF, leading_Poly.at(k).pow_CF )* //putting should have no effect if CF removed using leading
	//pow(2.0, leading_Poly.at(k).pow_2)* // changed 121218
	pow(TR, leading_Poly.at(k).pow_TR)*
	leading_Poly.at(k).int_part*
	leading_Poly.at(k).cnum_part;
	//cout << "cnum_num(Poly): Part k is "<< part_k << endl;
	res=res+part_k;
	}
	return res;
	}
	*/
	// To take the numerical value of a map
	//std::map< string, double > double_num(std::map< string, tr1::shared_ptr<Polynomial> > mem_map) {
	std::map< std::string, double > Col_functions::double_num( std::map< std::string, Polynomial > mem_map ) const{

	//cout << "double_num(std::map< string, tr1::shared_ptr<Polynomial> > mem_map) : Entering... " << endl;

	std::map< std::string, double > res;

	// Loop over entries
	//std::map< string, tr1::shared_ptr<Polynomial> >::iterator iter=mem_map.begin();
	std::map< std::string, Polynomial >::iterator iter=mem_map.begin();

	for( iter = mem_map.begin(); iter != mem_map.end(); ++iter) {
	//cout << "Hi " << endl;
	// Insert pair of string and the leading versions of the Polynomial
	double_num( (iter->second) );
	res.insert(std::make_pair( iter->first, double_num((iter->second)) ));

	}
	//cout << "double_num(std::map< string, tr1::shared_ptr<Polynomial> > mem_map) : Done." << endl;

	return res;
	}

	// Numerically evaluates a Polynomial for Nc=3,
	// and stores in the format of a Polynomial with only one term with only a numerical part
	Polynomial Col_functions::Polynomial_cnum_num( const Polynomial & Poly ) const{

	// Store content in numerical part of the Monomial
	Monomial Mon;
	Mon.cnum_part = cnum_num(Poly);

	Polynomial res;
	res = res * Mon;

	return res;
	}

	// Numerically evaluates a Poly_vec (vector of Polynomial) for Nc=3
	cvec Col_functions::cnum_num( const Poly_vec & Pv ) const{

	// To contain the numerical result
	cvec res;

	// Loop over Polynomials in the vector, and add the numerical value
	// to the vector to return
	for (uint p = 0; p < Pv.size(); p++) {
	res.push_back(cnum_num(Pv.at(p)));
	}
	return res;
	}

	// Numerically evaluates a Poly_vec
	dvec Col_functions::double_num( const Poly_vec & Pv ) const{
	//cout << "double_num(Pv): Incoming Pv \n" << Pv;
	// To contain the numerical result
	dvec res;

	// Loop over Polynomials in the vector, and add the numerical value
	// to the vector to return
	for (uint p = 0; p < Pv.size(); p++) {
	res.push_back(double_num( Pv.at(p)) );
	}
	return res;
	}

	// Numerically evaluates a Poly_matr
	dmatr Col_functions::double_num( const Poly_matr & Pm ) const{
	// To contain the numerical result
	dmatr res;

	// Loop over Poly_vecs in the vector, and add the numerical value
	// to the vector to return
	for (uint pv = 0; pv < Pm.size(); pv++) {
	res.push_back(double_num( Pm.at(pv).pv ) );
	}
	return res;
	}

	// Returns a double value. The argument is a vector of pointers to Polynomials.
	dvec Col_functions::double_num( const std::vector<std::tr1::shared_ptr<Polynomial> > & Pv ) const{
	//cout << "double_num(Pv): Incoming Pv \n" << Pv;
	// To contain the numerical result
	dvec res;

	// Loop over Polynomials in the vector, and add the numerical value
	// to the vector to return
	for (uint p = 0; p < Pv.size(); p++) {
	std::tr1::shared_ptr<Polynomial> the_pointer=Pv.at(p);
	// Want double of the polynomial which the pointer points at
	res.push_back( double_num(*the_pointer) );
	}
	return res;
	}


	// Numerically evaluates a Poly_vec (vector of Polynomial)
	// and stores in the form of a Polynomial
	Poly_vec Col_functions::Poly_vec_cnum_num( const Poly_vec & Pv) const{

	// To contain the result
	Poly_vec res;

	// Loop over Polynomials in the vector, put each polynomial to its numerical value
	for (uint p = 0; p < Pv.size(); p++) {
	res.push_back( Polynomial_cnum_num( Pv.at( p ) ) );
	}
	return res;
	}

	// Numerically evaluates a Poly_matr (vector of Poly_vec)
	// and stores in the form of a Poly_matr
	Poly_matr Col_functions::Poly_matr_cnum_num( const Poly_matr & Pm ) const {

	// To contain the result
	Poly_matr res_matr;

	// Loop over Polynomials in the matrix
	// and change each Polynomial to its numerical version
	for (uint v = 0; v < Pm.size(); v++) {
	res_matr.push_back( Poly_vec_cnum_num( Pm.at( v ).pv ));
	}

	return res_matr;
	}


	// Numerically evaluates a Poly_matr (vector of Poly_vec).
	cmatr Col_functions::cnum_num( const Poly_matr & Pm ) const{

	// To contain the numerical result
	cmatr res;

	// Loop over Polynomials in the vector, and add the numerical value
	// to the vector to return
	for (uint v = 0; v < Pm.size(); v++) {
	res.push_back(cnum_num( Pm.at(v).pv ));
	}
	return res;
	}



	dmatr Col_functions::double_num( const std::vector<std::vector<std::tr1::shared_ptr<Polynomial> > > & Pm ) const{
	//cout << " double_num(Pm): Incoming matrix \n" << Pm;

	// To contain the numerical result
	dmatr res;

	// Loop over Polynomials in the vector, and add the numerical value
	// to the vector to return
	for (uint v = 0; v < Pm.size(); v++) {
	res.push_back(double_num(Pm.at(v)));
	}
	return res;
	}


	/*
	// Numerically evaluates a poly_matr (vector of Poly_vec) for Nc=3
	cmatr leading_cnum_num(poly_matr Pm) {

	// To contain the numerical result
	cmatr res;
	res=cnum_num( leading( Pm ) );
	return res;
	}
	*/



	// A function to rename the indices in two Col_strs, such that in the first
	// they are called 1,2,3..., and in the second the relative order is kept
	std::pair<Col_str,Col_str> Col_functions::rename_indices( const Col_str & Cs1, const Col_str & Cs2 ) const{
	std::cout << "rename_indices: incoming Col_strs" << std::endl << Cs1 << Cs2;

	// To contain the new Cs1
	Col_str Cs1_new(Cs1);

	// To contain the new Cs2
	Col_str Cs2_new=(Cs2);

	// Loop over Quark_lines in Cs1
	int new_ind=0; // The new index
	for(uint i=0; i< Cs1_new.cs.size(); i++ ){
	//cout << i << endl;
	//cout.flush();
	// Loop over indices in the Quark_line
	for(uint j=0; j< Cs1_new.cs.at(i).ql.size(); j++ ){
	new_ind++;
	// The old index in Cs1, to be replaced also in Cs2
	int old_ind=Cs1_new.cs.at(i).ql.at(j);
	// Replace index in Cs1
	Cs1_new.cs.at(i).ql.at(j)=new_ind;
	// Locate and replace index in Cs2
	std::vector<int> place_in_2= Cs2.find_parton( old_ind );
	Cs2_new.cs.at( place_in_2.at(0) ).ql.at( place_in_2.at(1) )=new_ind;
	}
	}

	//cout << "rename_indices: outgoing Col_strs" << endl << Cs1 << Cs2_new;
	return std::make_pair(Cs1_new,Cs2_new);
	}


	// A function to rename the indices in two Col_amps, such that in the first
	// Col_str in the first Col_amp they are called 1,2,3...,
	// and in the other Col_strs the relative order is kept.
	std::pair<Col_amp,Col_amp> Col_functions::rename_indices( Col_amp Ca1, Col_amp Ca2 ) const{
	//cout << "rename_indices: incoming Col_strs" << endl << Cs1 << Cs2;

	// To contain the new Ca2
	Col_amp Ca2_new=Ca2;
	// The Col_str to define the ordering
	Col_str Cs1=Ca1.ca.at(0);
	// Loop over Quark_lines in Cs1
	int new_ind=0; // The new index
	for(uint i=0; i< Cs1.cs.size(); i++ ){
	//cout << i << endl;
	//cout.flush();
	// Loop over indices in the Quark_line
	for(uint j=0; j< Cs1.cs.at(i).ql.size(); j++ ){
	new_ind++;
	// The old index in Cs1, to be replaced also in Cs2
	int old_ind=Cs1.cs.at(i).ql.at(j);
	// Replace index in Cs1
	Cs1.cs.at(i).ql.at(j)=new_ind;
	// Locate and replace index in all other Col_strs
	// Do the replacement in all Col_strs in Ca1
	for(uint csi=0; csi < Ca1.ca.size(); csi++){
	std::vector<int> place= Ca1.ca.at(csi).find_parton( old_ind );
	Ca1.ca.at(csi).cs.at( place.at(0) ).ql.at( place.at(1) )=new_ind;

	}
	// Do the replacement in all Col_str in Ca2
	for(uint csi=0; csi < Ca2.ca.size(); csi++){
	std::vector<int> place= Ca2.ca.at(csi).find_parton( old_ind );
	Ca2_new.ca.at(csi).cs.at( place.at(0) ).ql.at( place.at(1) )=new_ind;
	}
	}
	}

	//cout << "rename_indices: outgoing Col_strs" << endl << Cs1 << Cs2_new;
	return std::make_pair(Ca1,Ca2_new);
	}


	// Function for calculating scalar products between Col_amps
	Polynomial Col_functions::scalar_product( const Col_amp & Ca1 , const Col_amp & Ca2 ) const{
	//cout <<"scalar_product(Ca1, Ca2): Entering " << endl;

	// To contain the result
	Col_amp Ca_res;

	// Contract the quarks
	Ca_res = contract_quarks(Ca1, Ca2);


	// Look for simple simplifications
	Ca_res.simplify();
	//Ca_res.simplify(); // why twice

	// Contract the gluons
	Ca_res = contract_all_gluons(Ca_res);
	//cout << "scalar_product: After contract_all_gluons(Cs_tmp ) " << Ca_res << endl;

	//Ca_res.remove_1_rings();
	//Ca_res.remove_0_rings();

	if (!Ca_res.ca.empty()) {
	// if( non_contracted ){
	std::cerr << "Col_functions::scalar_product: terminating due to non-contracted indices."
	<< std::endl;
	std::cerr << "The col_amp is " << Ca_res << std::endl;
	std::cerr.flush();
	throw -1;
	}

	else
	return Ca_res.Scalar;
	}


	// Alternative function for calculating scalar products between Col_amps
	Polynomial Col_functions::scalar_product2( const Col_amp & Ca1 , const Col_amp & Ca2 ) const{
	//cout <<"scalar_product(Ca1, Ca2): Entering " << endl;

	// To contain the result
	Polynomial Poly_res;
	Poly_res=Poly_res*0;

	// Loop over Col_strs
	for (uint cs1=0; cs1< Ca1.size(); cs1++){
	for (uint cs2=0; cs2< Ca2.size(); cs2++){

	// Calculate scalar product of Col_sts
	Poly_res = Poly_res+scalar_product( Ca1.at(cs1), Ca2.at(cs2) );
	}
	}
	double num_before=double_num(Poly_res);
	Poly_res.simplify();
	double num_after=double_num(Poly_res);

	if( abs(num_before - num_after) > accuracy ){
	std::cout <<"scalar_product2(Ca1, Ca2): Before simplification " << cnum_num(Poly_res) << std::endl;
	std::cout <<"scalar_product2(Ca1, Ca2): After simplification " << cnum_num(Poly_res) << std::endl;
	}
	return Poly_res;
	}


	// Function for calculating scalar product matrix
	dmatr Col_functions::scalar_product_matrix_mem( const Col_amp & Basis ) const{

	// To contain the resulting scalar product matrix
	std::vector<std::vector<std::tr1::shared_ptr<Polynomial> > > sp_m;

	// For remembering already calculated topologies
	std::map< std::string, std::tr1::shared_ptr<Polynomial> > mem_map;

	// Allocate memory as scalar product matrices can get large
	sp_m.reserve(Basis.ca.size());

	// Loop over basis vectors in Basis
	for ( uint i = 0; i < Basis.ca.size(); i++ ) {
	std::vector<std::tr1::shared_ptr<Polynomial> > rowi;

	//cout << "scalar_product_matrix_mem: calculating row " << i + 1 << " of "
	// << Basis.size() << endl;

	// Loop over basis vectors in Basis
	for ( uint j = 0; j < Basis.ca.size(); j++ ) {
	// To contain the ij-th entry
	std::tr1::shared_ptr<Polynomial> ijEntry;

	// Rename indices, and make string of new Col_strs, to use in map
	std::pair<Col_str, Col_str> Css_new = rename_indices(Basis.ca.at(i),
	Basis.ca.at(j));
	Col_str Cs1 = Css_new.first;
	Col_str Cs2 = Css_new.second;
	std::ostringstream Cs_string;
	Cs_string << Cs1 << Cs2;

	// Calculate element ij in scalar product matrix
	// If there are quarks involved this is straight forward
	if (Basis.ca.size() > 0 && !Basis.ca.at(0).gluons_only()) {
	// If this has scalar product has occurred before, reuse old value
	if ( mem_map.count(Cs_string.str()) > 0 ) {
	ijEntry = mem_map[Cs_string.str()];
	}
	// Otherwise, calculate value and save
	else {
	Polynomial p = scalar_product(Cs1, Cs2);
	ijEntry = shared_ptr<Polynomial> (new Polynomial(p));
	mem_map[Cs_string.str()] = ijEntry;
	}
	}
	// If there are only gluons involved the complex conjugate has to be
	// added (or subtracted) to consider the physical states
	// The sign is given by (-1)^Ng for Ng gluons, see arXiv:0906.1121 (my paper)
	if (Basis.ca.size() > 0 && Basis.ca.at(i).gluons_only()) {
	// The number of gluons in the basis stated
	// only correct and only used if gluons only
	uint Ng = Basis.ca.at(i).n_gluon();

	// The sign of the interference, (-1)^Ng
	int sign = (Ng % 2 ? -1 : 1);

	// We have 2 contributions, from \|amp\|^2 and \|conjugate(amp)\|^2 giving 2*Poly1,
	// The interferene comes later...
	std::tr1::shared_ptr<Polynomial> Poly1;
	Css_new = rename_indices(Basis.ca.at(i), Basis.ca.at(j));
	Cs1 = Css_new.first;
	Cs2 = Css_new.second;
	Cs_string.str().clear();
	Cs_string << Cs1 << Cs2;

	// Reuse value if calculated before
	if (mem_map.count(Cs_string.str()) > 0) {
	Poly1 = mem_map[Cs_string.str()];
	}
	// Otherwise, calculate value and save
	else {
	Polynomial p = scalar_product(Cs1, Cs2);
	Poly1 = shared_ptr<Polynomial> (new Polynomial(p));
	mem_map[Cs_string.str()] = Poly1;
	}

	// For the second part we need the conjugate
	// Rename indices, and make string of new Col_strs
	Col_str Basis_j_conj=Basis.ca.at(j);
	Basis_j_conj.conjugate();
	std::pair< Col_str, Col_str > Css_new = rename_indices(Basis.ca.at(i), Basis_j_conj );
	Cs1 = Css_new.first;
	Cs2 = Css_new.second;
	Cs_string.str().clear();
	Cs_string << Cs1 << Cs2;

	// For the interference
	std::tr1::shared_ptr<Polynomial> Poly2;
	// Reuse value if calculated before
	if (mem_map.count(Cs_string.str()) > 0) {
	//cout << "Already calculated sp!" << endl;
	Poly2 = mem_map[Cs_string.str()];
	}
	// Otherwise, calculate value and save
	else {
	Polynomial p = scalar_product( Basis.ca.at(i), Basis_j_conj );
	Poly2 = shared_ptr<Polynomial> (new Polynomial(p));
	mem_map[Cs_string.str()] = Poly2;
	}
	Polynomial p = 2 * (Poly1) + sign 2 * (*Poly2);
	ijEntry = shared_ptr<Polynomial> (new Polynomial(p));
	(*ijEntry).simplify();
	}
	rowi.push_back(ijEntry);
	}
	sp_m.push_back(rowi);
	}


	// Leading terms version
	dmatr diag_matr= leading(sp_m);
	// Verifying that the leading terms only sit on the diagonal
	bool diagonal=check_diagonal(diag_matr);
	if( !diagonal ) {throw -1;
	std::cerr <<"scalar_product_matrix_mem: leading terms are not diagonal" << std::endl;
	std::cerr.flush();
	}
	write_out_spm( diag_matr, Basis , true );
	// Clear unneeded matrix as it may be large
	diag_matr.clear();

	// Numerical version for checking symmetry and saving
	dmatr matr= double_num(sp_m);

	// Verifying that the matrix is symmetric to accuracy
	bool sym=check_symmetry( matr );
	if( !sym ) {
	std::cerr <<"scalar_product_matrix_mem: scalar product matrix not symmetric" << std::endl;
	std::cerr.flush();
	throw -1;
	}
	// Also save the numerical versions to files for future usage
	// Full Nc version
	write_out_spm( matr, Basis , false );


	return matr;
	}

	// Function for calculating scalar product matrix
	dmatr Col_functions::scalar_product_matrix_mem2( const Col_amp & Basis, const int n_loop ) const{

	// To contain the resulting scalar product matrix
	dmatr matr;
	dmatr matr_l;

	// For remembering already calculated topologies
	std::map< std::string, Polynomial> mem_map;

	std::cout << "scalar_product_matrix_mem2: calculating topologies... " << std::endl;
	// Loop over basis vectors in Basis once, this should give all topologies unless gluos only
	// and loop level

	if( n_loop >0 or Basis.n_quark()>0 ) {
	std::cout << "scalar_product_matrix_mem2: Loop level calculation or >1 q qbar, calling scalar_product_matrix_mem instead." << std::endl;
	return scalar_product_matrix_mem( Basis);
	}
	// Loop over basis vectors in Basis once, this should give all topologies unless gluons only
	// and loop level
	uint max = 1;
	for (uint i = 0; i <max; i++) {

	for (uint j = 0; j < Basis.size(); j++) {
	// Rename indices, and make string of new Col_strs, to use in map
	std::pair< Col_str, Col_str > Css_new = rename_indices(Basis.ca.at(i),
	Basis.ca.at(j));
	Col_str Cs1 = Css_new.first;
	Col_str Cs2 = Css_new.second;
	std::ostringstream Cs_string;
	Cs_string << Cs1 << Cs2;

	// Add result to map
	//if (!(mem_map.count(Cs_string.str()) > 0)) {
	if ((mem_map.find(Cs_string.str()) == mem_map.end())) {
	//cout << "adding term for" << endl << Cs_string.str();
	Polynomial p; // to contain result

	// In the case of both quarks and gluons
	if (!Basis.gluons_only()) {
	p = scalar_product(Cs1, Cs2);
	}
	// For gluons with implicitly added topology
	// In the case of only gluons the conjugated topology is needed as well
	// Note that the full result is saved in the local map
	else {
	uint Ng = Basis.n_gluon();
	int sign = (Ng % 2 ? -1 : 1);
	Col_str Cs2_conj=Cs2;
	Cs2_conj.conjugate();
	p = 2 * scalar_product(Cs1, Cs2) + 2 * sign * scalar_product(Cs1,
	Cs2_conj);
	p.simplify();

	}
	if( i>0 )
	std::cout << "the combination " << Cs_string.str() << " was not found before" << std::endl;
	mem_map.insert(make_pair(Cs_string.str(), p));

	//cout << "The memmap at place " << mem_map[Cs_string.str()] ;
	}
	// For gluons, also the conjugated topology can appear, add this as well to map
	if (Basis.gluons_only()) {

	//Cs2=conjugate(Basis.ca.at(j));
	Cs2=Basis.ca.at(j);
	Cs2.conjugate();
	Cs2.simplify();
	Css_new = rename_indices( Basis.ca.at(i), Cs2 );
	//cout << "the conjugate(Basis.ca.at(j)) " << Cs2 << endl;

	Cs1 = Css_new.first;
	Cs2 = Css_new.second;
	Col_str Cs2_conj=Cs2;
	Cs2_conj.conjugate();
	std::ostringstream Cs_string;
	Cs_string << Cs1 << Cs2;
	if (!(mem_map.count(Cs_string.str()) > 0)) {
	Polynomial p; // to contain result
	uint Ng = Basis.n_gluon();
	int sign = (Ng % 2 ? -1 : 1);
	p = 2 * scalar_product(Cs1, Cs2) + 2 * sign * scalar_product(Cs1, Cs2_conj);
	p.simplify();
	mem_map.insert(make_pair(Cs_string.str(), p));
	}
	}

	} // End of j-loop
	} // End of i-loop
	std::cout << "done. ";
	std::cout.flush();
	std::cout << mem_map.size() << " possible topologies. " << std::endl;

	// Now that the map is found construct numerical map, and the leading Nc map
	std::cout << "scalar_product_matrix_mem2: Finding numerical entries " << std::endl;
	std::map< std::string, double> num_mem_map = double_num(mem_map);
	std::cout << "done " << std::endl;
	std::cout << "scalar_product_matrix_mem2: Finding leading entries " << std::endl;
	std::map<std::string, double> num_leading_mem_map = double_num(leading(mem_map));
	std::cout << "done " << std::endl;
	std::cout << "There was " << num_leading_mem_map.size()
	<< " numerical leading topologies. " << std::endl;

	// Calculating numerical scalar product matrices
	for (uint i = 0; i < Basis.ca.size(); i++) {

	// To contain vectors in matrices
	//vector< Polynomial > rowi;
	dvec rowi;
	dvec rowi_l;

	std::cout << "scalar_product_matrix_mem2: calculating rows " << i + 1 << " of "
	<< Basis.size() << std::endl;

	// Loop over basis vectors in Basis
	for (uint j = 0; j < Basis.ca.size(); j++) {
	// To contain the ij-th entry
	double ijEntry;
	double ijEntry_l;

	// Rename indices, and make string of new Col_strs, to use in map
	std::pair<Col_str, Col_str> Css_new = rename_indices(Basis.ca.at(i),
	Basis.ca.at(j));
	Col_str Cs1 = Css_new.first;
	Col_str Cs2 = Css_new.second;
	std::ostringstream Cs_string;
	Cs_string << Cs1 << Cs2;

	// Add to ijth entries, in case of both quarks and gluons
	if ((mem_map.count(Cs_string.str()) > 0)) {
	ijEntry = num_mem_map[Cs_string.str()];
	ijEntry_l = num_leading_mem_map[Cs_string.str()];
	//cout << "top found" << endl;
	} else {
	std::cerr
	<< "scalar_product_matrix_mem2: All scalar product topologies should have been found by now."
	<< std::endl;
	std::cerr.flush();
	throw -1;
	}

	rowi.push_back(ijEntry);
	rowi_l.push_back(ijEntry_l);
	} // End of j loop

	matr.push_back(rowi);
	matr_l.push_back(rowi_l);
	} // End of i loop

	// Verifying that the matrix is symmetric to accuracy
	bool sym = check_symmetry(matr);
	if (!sym)
	throw -1;
	// Also save the numerical versions to files for future usage
	std::cout << "scalar_product_matrix_mem2: Writing out scalar product matrix... "
	<< std::endl;
	write_out_spm(matr, Basis, false);

	// Verifying that the leading terms only sit on the diagonal
	bool diagonal = check_diagonal(matr_l);
	if (!diagonal)
	throw -1;
	std::cout
	<< "scalar_product_matrix_mem2: Writing out leading contribution to scalar product matrix... "
	<< std::endl;
	write_out_spm(matr_l, Basis, true);

	return matr;
	}


	Polynomial Col_functions::scalar_product( const Col_str & Cs1, const Col_str & Cs2 ) const{

	//cout << "scalar_product(Cs1, Cs2): has the arguments " << Cs1 << " " << Cs2 << endl;

	Col_str Cs_tmp;
	Col_amp Ca_tmp;


	// Contract the quarks
	Cs_tmp=contract_quarks( Cs1, Cs2 );
	//cout << "scalar_product: After contraction of quarks" << Cs_tmp << endl;
	std::cout.flush();

	// Contract the gluons
	Ca_tmp=contract_all_gluons(Cs_tmp );
	//cout << "scalar_product: After contract_all_gluons(Cs_tmp ) " << Ca_tmp << endl;


	//if ( Ca_tmp.ca.empty() ) cout << "empty ca" << endl;
	//cout << "Ca_tmp.poly "<< Ca_tmp.poly <<endl;
	//cout << "scalar_product(Cs1,Cs2): To return, " << Ca_tmp.Poly << endl;
	if ( !Ca_tmp.ca.empty() ){
	// if( non_contracted ){
	std::cerr << "scalar_product: terminating due to non-contracted indices." <<std::endl;
	std::cerr << "The col_amp is " << Ca_tmp << std::endl;
	throw -1;
	}


	else return Ca_tmp.Scalar;
	}

	// Function to test if a numerical (scalar product matrix) is symmetric
	bool Col_functions::check_symmetry( const dmatr & matr ) const{
	std::cout << "check_symmetry( matr ): Numerically verifying symmetry of matrix...";
	if( matr.empty() ==0 ){
	std::cout << "The numerical matrix is empty. " << std::endl;
	return true;
	}

	// Verifying that the matrix is symmetric to accuracy "accuracy",
	// if not sym is put to false
	bool sym=true;
	// Loop over basis vectors in Basis
	for (uint i = 0; i < matr.size(); i++) {
	std::vector<Polynomial> rowi;
	// Loop over basis vectors in Basis
	for (uint j = 0; j <=i; j++) {
	if ((abs(matr.at(i).at(j) / matr.at(j).at(i) - 1.0) > accuracy) && (matr.at(i).at(j) > accuracy) && (matr.at(j).at(i) > accuracy)) {
	sym=false;
	std::cerr
	<< "check_symmetry( matr ): Error, the resulting scalar product matrix is not symmetric. \n "
	<< "Element " << i << "," << j << ": "
	<< matr.at(i).at(j) << ", Element " << j << "," << i
	<< ": " << matr.at(j).at(i) << std::endl
	<< "This indicates an error in calculation of scalar products. " << std::endl;
	}
	}
	}
	std::cout << " done." << std::endl;
	return sym;
	}

	// Function to test if a numerical (scalar product matrix) is diagonal
	bool Col_functions::check_diagonal(const dmatr & matr) const{
	std::cout << "check_diagonal( matr ): Numerically checking if the matrix is diagonal...";

	// Verifying that the matrix is diagonal
	// if not sym is put to false
	bool diag=true;
	// Loop over basis vectors in Basis
	for (uint i = 0; i < matr.size(); i++) {
	std::vector<Polynomial> rowi;
	// Loop over basis vectors in Basis
	for (uint j = 0; j <=i; j++) {
	if (matr.at(i).at(j)!=0.0 && i!=j) {
	diag=false;
	std::cerr
	<< "check_diagonal( matr ): Warning, the matrix is not diagonal. \n "
	<< "Element " << i << "," << j << ": "
	<< matr.at(i).at(j) << std::endl
	<< "This indicates an error in calculation of scalar products. " << std::endl;

	}
	}
	}
	std::cout << " done." << std::endl;
	return diag;
	}

	/*
	// A function for calculating the scalar products between the basis vectors
	// (stored in a Col_amp)
	Poly_matr Col_functions::scalar_product_matrix( const Col_amp & Basis ) const{

	// To contain the resulting scalar product matrix
	Poly_matr sp_m;

	// Loop over basis vectors in Basis
	for( uint i=0; i< Basis.ca.size(); i++){
	std::vector< Polynomial > rowi;

	// Loop over basis vectors in Basis
	for( uint j=0; j< Basis.ca.size(); j++){
	//cout << "scalar_product_matrix(Col_amp):: calculating for i " << i << " j " << j << endl;
	//cout.flush();
	// Calculate element ij in scalar product matrix
	// If there are quarks involved this is straight forward
	Polynomial ijEntry;
	if( Basis.ca.size()>0 && !Basis.ca.at(i).gluons_only() ){
	ijEntry=scalar_product(Basis.ca.at(i), Basis.ca.at(j));
	}
	// If there are gluons involved the complex conjugate has to be
	// added or subtracted to consider the physical states as
	// the conjugate, "anti-cyclic", amplitude is implicit
	// The sign is given by (-1)^Ng for Ng gluons, see arXiv:0906.1121 (my paper)
	if( Basis.ca.size()>0 && Basis.ca.at(i).gluons_only() ){
	// The number of gluons in the basis stated
	// only correct and only used if gluons only
	uint Ng=Basis.ca.at(i).n_gluon();

	// The sign of the interference, (-1)^Ng
	int sign=(Ng % 2 ? -1:1);
	//cout << "The sign is " << sign << " for Ng= " << Ng <<endl;
	// If only gluons the basis states are composed of 2 color structures
	// Add together to get scalar product matrix for basis
	// Using the realness of the vector space
	ijEntry=2*scalar_product( Basis.ca.at(i), Basis.ca.at(j) )
	+sign2scalar_product( Basis.ca.at(j), conjugate( Basis.ca.at(i) ) );
	ijEntry.simplify();
	}
	rowi.push_back( ijEntry );
	}
	sp_m.push_back(rowi);
	}

	// For saving and checking symmetry of numerical version
	dmatr matr= double_num(sp_m);

	// Verifying that the matrix is symmetric to accuracy
	bool sym=check_symmetry( matr );
	if( !sym ) throw -1;

	// Also save the numerical versions to files for future usage
	// Full Nc version
	write_out_spm( matr, Basis , false );

	// Leading terms version
	matr= double_num( leading(sp_m) );
	write_out_spm( matr, Basis , true );

	// Verifying that the leading terms only sit on the diagonal
	bool diagonal=check_diagonal(matr);
	if(!diagonal) throw -1;

	return sp_m;
	}
	*/
	// Functions for generating the n_partons random amplitudes in [-1/2,1/2]
	// for emission from n_partons.
	cvec Col_functions::generate_random_amplitudes( int n_p ) const{

	cvec amps;

	// Loop over partons which can radiate
	for(int i=0; i< n_p; i++ ){
	// Random number in [-0.5,0.5]
	double r=1.0;
	r=rand()/float(RAND_MAX)-0.5;
	amps.push_back(r);
	}
	return amps;
	}


	Poly_vec Col_functions::decompose( const Col_amp & Ca, const Col_amp & Basis ) const{
	//cout << "decompose: Incoming Basis and Col_amp " << Basis << Ca << endl;
	// To contain the decomposed vector
	Poly_vec Decv;
	Polynomial Zero;
	Zero=Zero*0;

	// Initially set all components of the vectors to 0
	for (uint m2 = 0; m2 < Basis.ca.size(); m2++){
	Decv.push_back(Zero);
	}
	//cout<<"decompose: initial Decv "<< Decv << "with size "<< Decv.size() << endl;

	// Loop over Cs in Ca and check which basis vector they equal
	for (uint m1 = 0; m1 < Ca.ca.size(); m1++) {
	for (uint m2 = 0; m2 < Basis.ca.size(); m2++) {
	bool found=false;
	if (Ca.ca.at(m1).cs == Basis.ca.at(m2).cs) {
	found=true;
	//cout << m2 << " Parts " << Decv.at(m2) << " and " << Ca.ca.at(m1).Poly << endl;
	//if( Ca.ca.at(m1).Poly.empty() ) cout << m2 << "Poly was empty" << endl;
	Decv.at(m2)=Decv.at(m2)+Ca.ca.at(m1).Poly;
	//cout.flush();
	//if(found) cout << m2 << " result " << Decv.at(m2) << " num " << double_num(Decv.at(m2)) << endl;
	//cout.flush();

	}

	if(found) break; // The tensor is ONE of the tensors in the Basis
	}
	}
	//cout << "decompose: final Decv "<< Decv << "with size "<< Decv.size() << endl;
	return Decv;
	}


	// Function for calculating scalar products
	// given the information about the basis and the scalar product matrix in the basis
	Polynomial Col_functions::scalar_product( const Col_amp & Ca1, const Col_amp & Ca2, const Col_amp & Basis, const Poly_matr & sp_matr ) const{

	// To contain the resulting Polynomial
	Polynomial Poly_res;
	Poly_res=Poly_res*0;

	// Decompose the Col_amps
	Poly_vec Polyv1=decompose(Ca1, Basis);
	Polyv1.conjugate();// conjugation added 130530
	Poly_vec Polyv2=decompose(Ca2, Basis);

	//std::cout << "scalar_product(Ca1,Ca2,Basis,sp_matr): Decomposition done."<< std::endl;
	// Then add contributions
	for (uint m1=0; m1< Basis.ca.size(); m1++){
	// Diagonal terms
	Poly_res=Poly_res+Polyv1.at(m1) Polyv2.at(m1) sp_matr.at(m1).at(m1);
	for (uint m2=0; m2< m1; m2++){
	// Other terms, use symmetry of scalar product matrix
	//Poly_res=Poly_res+ 2Polyv1.at(m1) Polyv2.at(m2) *sp_matr.at(m1).at(m2); how could this be right
	Poly_res=Poly_res+ ( Polyv1.at(m1) Polyv2.at(m2)+Polyv1.at(m2) Polyv2.at(m1) ) *sp_matr.at(m1).at(m2);

	}
	}
	return Poly_res;
	}

	// Function for calculating scalar products numerically
	// given the information about the basis and the scalar product matrix in the basis
	Polynomial Col_functions::scalar_product_cnum_num( const Col_amp & Ca1, const Col_amp & Ca2, const Col_amp & Basis, const Poly_matr & sp_matr ) const{

	// To contain the resulting Polynomial
	Polynomial Poly_res;
	Poly_res=Poly_res*0;

	// Decompose the Col_amps
	Poly_vec Polyv1=decompose(Ca1, Basis);
	// added complex conjugation 130525
	Polyv1.conjugate();

	Poly_vec Polyv2=decompose(Ca2, Basis);
	//std::cout << "scalar_product(Ca1,Ca2,Basis,sp_matr): Decomposition done."<< std::endl;
	// Then add contributions
	for (uint m1=0; m1< Basis.ca.size(); m1++){
	// Diagonal terms
	Poly_res=Poly_res+ Polyv1.at(m1) Polyv2.at(m1) sp_matr.at(m1).at(m1);
	for (uint m2=0; m2< m1; m2++){
	// Other terms, use symmetry of scalar product matrix
	//Poly_res=Poly_res+ 2Polyv1.at(m1) Polyv2.at(m2) *sp_matr.at(m1).at(m2); // How could this be right?
	// added complex conjugation 130525
	Poly_res+= ( Polyv1.at(m1) Polyv2.at(m2) +Polyv1.at(m2) Polyv2.at(m1) ) *sp_matr.at(m1).at(m2);

	}
	}

	return Poly_res;

	}


	// Function for calculating scalar products
	// given the information about the basis in the scalar product matrix in the basis
	// This is an attempt to speed up
	cnum Col_functions::scalar_product_cnum_num( const Col_amp & Ca1, const Col_amp & Ca2, const Col_amp & Basis, const dmatr & sp_matr ) const{

	// To contain the resulting Polynomial
	cnum res=0;
	uint basis_size=Basis.ca.size();
	//time_t t1, t2, t3;
	//t1=clock( );
	// Decompose the Col_amps
	Poly_vec Polyv1, Polyv2;
	Polyv1=decompose(Ca1, Basis);
	if( &Ca1==&Ca2 ) Polyv2=Polyv1;
	else Polyv2=decompose(Ca2, Basis);
	//cout << "scalar_product_cnum_num(Ca1, Ca2, Basis, sp_matr)" << "\n Polyv1: " << Polyv1 << endl;

	// Conjugate first arg, added 130525
	Polyv1.conjugate();


	// Make ordinary vectors to speed up multiplication
	cvec v1, v2;
	v1.reserve( basis_size );
	for (uint i=0; i< basis_size; i++){
	v1.push_back( cnum_num( Polyv1.at(i) ) );
	};
	//if( &Ca1==&Ca2 ) v2=v1;
	if (false) v2=v1;
	else{
	v2.reserve(basis_size);
	for (uint i=0; i< basis_size; i++) v2.push_back( cnum_num( Polyv2.at(i) ) );
	}
	//cout << "scalar_product_cnum_num(Ca1, Ca2, Basis, sp_matr) " << "\n v1: " << v1 << endl;
	//t2=clock();
	//cout << "scalar_product_cnum_num(Ca1,Ca2,Basis, cmatr): Decomposition done at time: " << t2-t1 << endl;

	// Then add contributions
	cnum v1m1, v2m1;
	dvec row;
	for (uint m1=0; m1< basis_size; m1++){
	v1m1=v1.at(m1);
	v2m1=v2.at(m1);
	row=sp_matr.at( m1 );
	// Diagonal terms
	res=res+ v1m1 v2.at(m1) row.at(m1);
	for (uint m2=0; m2< m1; m2++){
	//res=res+ v1.at(m1) v2.at(m2) sp_matr.at(m1).at(m2);
	//res=res+ 2.0ind1v2.at(m2) *row.at(m2);// how could this be true
	res += (v1m1v2.at(m2) +v1.at(m2)v2m1 ) *row.at(m2); // changed to this 130524

	}
	}
	//t3=clock();
	//cout << "scalar_product_cnum_num(Ca1,Ca2,Basis, cmatr): Matrix multiplication done at time: " << t3-t2 << endl;
	return res;

	}


	// Connect n_q q and qbars in all n_q! ways, returns a basis (in the form of a Col_amp)
	// for the case of n_q q qbar pairs.
	Col_amp Col_functions::connect_quarks( int n_q ) const{
	Col_amp Basis;
	// If n_q=1, there is only one state
	if( n_q==1 ) {
	Col_str Cs_tmp("[{1,2}]");
	Basis=Basis+Cs_tmp;
	}
	else{ // For more than one q qbar pair, build up states iteratively
	// Start from the state with one less q qbar state
	Col_amp Old_bas=connect_quarks( n_q-1 );

	// Loop over old basis vectors
	for (uint old_v = 0; old_v < Old_bas.ca.size(); old_v++) {

	// Some, 1*(n_q-1)!, new states are obtained by just adding q_new qbar_new to the old states
	// Construct the new Ql with i.e. {q_new, qbar_new}
	Quark_line Ql_new;
	Ql_new.open = true;
	Ql_new.ql.push_back(2 * n_q - 1);
	Ql_new.ql.push_back(2 * n_q);
	Col_str New_state;
	New_state.cs.push_back(Ql_new);
	New_state.append( Old_bas.ca.at(old_v).cs );
	Basis=Basis+New_state;
	}

	// The other (n_q-1)*(n_q-1)! states are obtained by combining the new qbar with any old quark
	// Loop over old basis states
	for (uint old_v = 0; old_v < Old_bas.ca.size(); old_v++) {
	// Loop over old qbar indices
	for (int qbar_old = 2; qbar_old < 2 * n_q; qbar_old += 2) {
	// To contain the new basis vector
	Col_str New_state = Old_bas.ca.at(old_v);

	// Replace the index qbar_old with 2*nq
	New_state.replace(qbar_old, 2 * n_q);

	// Create the Quark_line {q_new, qbar_old}
	Quark_line Ql_new;
	Ql_new.open = true;
	Ql_new.ql.push_back(2 * n_q - 1);
	Ql_new.ql.push_back(qbar_old);

	// Add the Col_str from {q_new, qbar_old} to the new basis vector
	New_state.cs.push_back(Ql_new);

	// Add the new basis tensor to the basis
	Basis = Basis + New_state;
	}
	}
	}

	return Basis;
	}

	// Function for calculating scalar products
	// This function assumes that the Col_amps are normal ordered and that the full polynomial is
	// contained in the numerical factor
	// It also assumes that the BASIS in which sp_matr is calculated is the SAME as that of the Col_amp
	// and that the Col_amp is a linear combination of ALL states
	// The coefficients of the basis vectors can then simply be read off in the Polynomials multiplying the Col_strs
	// This is an attempt to speed up
	cnum Col_functions::scalar_product_num_fast( const Col_amp & Ca1, const Col_amp & Ca2, const dmatr & sp_matr ) const{

	// To contain the result
	cnum res=0;
	uint basis_size=Ca1.ca.size(); // Both Ca's should have same length
	//time_t t1, t2, t3;
	//t1=clock( );

	// Make ordinary vectors of numerical information in Polynomials, to speed up multiplication
	cvec v1, v2;
	v1.reserve(basis_size);
	for (uint i=0; i< basis_size; i++){
	// added conjugation of first arg 130524
	v1.push_back( conj( cnum_num( Ca1.ca.at(i).Poly ) ) );
	};
	if( &Ca1==&Ca2 ) v2=v1;
	else{
	v2.reserve(basis_size);
	for (uint i=0; i< basis_size; i++) v2.push_back( cnum_num( Ca1.ca.at(i).Poly ) );
	}

	//t2=clock();

	// Then add contributions
	cnum ind1;
	dvec row;
	for (uint m1=0; m1< basis_size; m1++){
	ind1=v1.at(m1);
	row=sp_matr.at(m1);
	// Diagonal parts
	res=res+ ind1 v2.at(m1) row.at(m1);
	for (uint m2=0; m2< m1; m2++){
	//res=res+ v1.at(m1) v2.at(m2) sp_matr.at(m1).at(m2);
	// Non diagonal parts
	res=res+ 2.0ind1 v2.at(m2) *row.at(m2);
	}
	}
	//t3=clock();

	return res;

	}

	// Calculates the scalar product between numerical (complex) amplitudes v1, v2 given the scalar product matrix
	// Warning: ordering of basis tensors must be the same
	cnum Col_functions::scalar_product_num_fast( const cvec & v1, const cvec & v2, const dmatr & spm_matr ) const {

	if(v1.size()!= v2.size()){
	std::cerr << "Col_functions::scalar_product_num_fast: Size of first vector "
	<< v1.size() << " does not agree with size of second vector "
	<< v2.size() << std::endl;
	throw v2;
	}
	if(v1.size()!= spm_matr.size()){
	std::cerr << "Col_functions::scalar_product_num_fast: Size of vectors "
	<< v1.size() << " do not agree with size of matrix "
	<< spm_matr.size() << std::endl;
	throw spm_matr;
	}
	uint basis_size=v1.size();

	// To contain the result
	cnum res=0;
	cnum tmp=0, tmp2=0;

	//time_t t1, t2;
	//t1=clock( );

	// Add contributions
	cnum v1m1, v2m1;
	for (uint m1=0; m1< basis_size; m1++){
	v1m1=v1.at(m1);
	v2m1=v2.at(m1);
	const dvec &row = spm_matr.at(m1);
	// Diagonal parts
	res += v1m1 v2.at(m1) row.at(m1);
	tmp = 0;
	tmp2=0;
	for (uint m2=0; m2< m1; m2++){
	//res=res+ v1.at(m1) v2.at(m2) sp_matr.at(m1).at(m2);
	//res=res+ 2.0ind1 v2.at(m2) *row.at(m2);
	//tmp += v2[m2] *row[m2]; // 40% quicker
	tmp += (v1m1v2[m2]+v2m1v1[m2]) *row[m2];

	//alt
	//tmp+=v2[m2]*row[m2];
	//tmp2+=v1[m2]*row[m2];


	}
	//res += 2.0v1m1tmp;
	res+=tmp;
	// alt
	//res+= v1m1tmp+v2m1tmp2;
	}
	//t2=clock();
	//cout << "scalar_product_num_fast(v1, v2, spm_matr): scalar prod calculated at time "
	// << t2-t1 << endl;
	return res;
	}

	// Calculates the scalar product between numerical (complex) amplitudes v1, v2 given the scalar product matrix
	// Warning: ordering of basis tensors must be the same
	// Function for calculating scalar products, assuming that only diagonal parts contribute
	// This is useful for calculating the leading Nc result
	// This function assumes that the Col_amps are normal ordered and that the full polynomial is
	// contained in the numerical factor
	// It also assumes that the BASIS in which sp_matr is calculated is the SAME as that of the Col_amp
	// and that the Col_amp is a linear combination of ALL states
	// The coefficients of the basis vectors can then simply be read off in the Polynomials multiplying the Col_strs
	// This is an attempt to speed up
	// Calculates the scalar product between numerical (complex) amplitudes v1, v2 given the scalar product matrix
	// Assumes that there are only diagonal contributions. This is useful for calculations in leading Nc limit
	// Warning: ordering of basis tensors must be the same
	cnum Col_functions::scalar_product_num_diagonal( const cvec & v1, const cvec & v2, const dmatr & spm_matr ) const {

	// To contain the result
	cnum res=0;

	// Check dimension
	uint basis_size=v1.size();

	//time_t t1, t2;
	//t1=clock( );

	// Then add contributions
	for (uint m1=0; m1< basis_size; m1++){
	// added conjugation of first arg 130524
	res=res+ conj( v1.at(m1) ) v2.at(m1) spm_matr.at(m1).at(m1);
	}
	//t2=clock();
	//cout << "scalar_product_num_diagonal(v1, v2, spm_matr): scalar prod calculated at time "
	// << t2-t1 << endl;
	return res;
	}

	// Compute the new basis states coming from one old vector when one gluon, g_new, is added.
	Col_amp Col_functions::add_one_gluon( const Col_str & Cs, int g_new, int n_loop ) const{
	//cout << "add_one_gluon( Col_str, int): Entering with arg" << Old_tens
	// << endl;

	// For storing the new basis
	Col_amp New_tensors;

	// Add the new gluon in all possible ways to the old Color structure
	// Loop over the Quark_lines
	for (uint ql = 0; ql < Cs.cs.size(); ql++) {
	// The old Quark_line, before insertion of the new gluon index
	Quark_line Old_Ql = Cs.cs.at(ql);
	Col_str New_tensor = Cs;

	// Special case of insertion of a gluon in a 2-ring in a gluons only basis
	// in this case only "half" the basis states for rings with >= 3 gluons are
	// generated. The rest are obtained by taking the indices in anti-cyclic order
	// i.e. by complex conjugating
	if ((Old_Ql.ql.size() == 2 && Cs.gluons_only())) {
	// Insert the new gluon after the existing gluons
	Quark_line New_Ql = Old_Ql;
	New_Ql.ql.push_back(g_new);
	// Replace the old Col_str with the new and add to the new basis states
	New_tensor.cs.at(ql) = New_Ql;
	New_tensors = New_tensors + New_tensor;
	} else { // ordinary case
	// Loop over (potential) insertion places in the Quark_lines, starting from the end
	for (int j = Old_Ql.ql.size(); j > 0; j--) {
	//cout << "add_one_gluon( Col_str, int): ql " << ql << " j " << j
	// << endl;
	// Special treatment of last place, insert here only if the ring is open
	// (the gluon index cannot take the place of the a quark index)
	if (Old_Ql.open && j == static_cast<int>(Old_Ql.ql.size()) ) j--;

	Quark_line New_Ql = Old_Ql;
	quark_line::iterator it = New_Ql.ql.begin() + j;
	New_Ql.ql.insert(it, g_new);

	//// Inserting element at place given by iterator itj
	//quark_line::iterator itj = cs.at(i).ql.begin() + j;
	//cs.at(i).ql.insert(itj, part_num);

	// Replace the old Col_str with the new and add to the new basis states
	New_tensor.cs.at(ql) = New_Ql;
	New_tensors = New_tensors + New_tensor;
	}
	}
	}
	//cout << "add_one_gluon( Col_str, int): to return " << New_tensors << endl;
	return New_tensors;
	}

	// Compute the basis states coming from the whole Col_amp if one gluon is added
	// If n_loop==0, only tree level states are constructed.
	Col_amp Col_functions::add_one_gluon( const Col_amp & Old_basis, int n_q, int n_g, int g_new, int n_loop ) const{
	//cout << "add_one_gluon( Col_amp, int): Entering with arg " << Old_bas
	// << endl;
	// For storing the new basis
	Col_amp New_bas;

	//cout << "add_one_gluon( Col_amp, int): current state " << New_bas << endl;

	// Add the new gluon to each of the previous color tensors
	for (uint t = 0; t < Old_basis.ca.size(); t++) {
	New_bas = New_bas + add_one_gluon(Old_basis.ca.at(t), g_new, n_loop);
	//cout << "add_one_gluon( Col_amp, int): current state " << New_bas << endl;
	}

	//cout << "add_one_gluon( Col_amp, int): Gluons added to old states done with result " << New_bas
	// << endl;

	// Create new states with 2-rings formed from taking the new gluon and any old gluon
	// these states can be created by taking the n_g-2 basis but use the indices
	// 1...i_ex-1, i_ex+1 n_g-1
	// for the gluons, that is exclude one index (to be combined with the index n_g in a 2-ring)
	// First, create the n_g-2 basis here, must have at least 4 gluons
	// Do this only is not tree-level mode, i.e. n_loop>0
	if ( n_loop>0 && (( n_q>0 && g_new-2*n_q>=2 ) or (n_q==0 && g_new >= 4) )) {
	Col_amp Bas_n_g_minus_2 = create_trace_basis(n_q, g_new-n_q*2-2, n_loop);

	// cout << "add_one_gluon( Col_amp, int): Bas_n_g_minus_2 " << Bas_n_g_minus_2
	// << endl;

	// Loop over indices to group with the new gluon index
	//for(int g_ex=2*n_q+1; g_ex <= n_g-2; g_ex++){
	for (int g_ex = 2 * n_q + 1; (g_ex <= g_new - 1 ); g_ex++) {
	//cout
	// << "add_one_gluon( Col_amp, int): Looping over excluded index, g_ex "
	// << g_ex << endl;
	// Create the 2-ring with g_ex and g_new
	Quark_line Ql2_ring;
	Ql2_ring.open = false;
	Ql2_ring.append(g_ex);
	Ql2_ring.append(g_new);
	//cout << "add_one_gluon( Col_amp, int): Created 2-ring " << Ql2_ring
	// << endl;

	// Loop over Col_str's=basis vectors
	for (uint bv = 0; bv < Bas_n_g_minus_2.ca.size(); bv++) {
	//cout << "add_one_gluon( Col_amp, int): Looping over bv " << bv
	// << endl;
	Col_str Cs_new; // To contain the new basis vecto
	Cs_new = Bas_n_g_minus_2.ca.at(bv);
	// Loop over the Quark_lines in the basis vectors
	for (uint ql = 0; ql < Bas_n_g_minus_2.ca.at(bv).cs.size(); ql++) {
	//cout << "add_one_gluon( Col_amp, int): Looping over ql "
	// << ql << endl;
	// Loop over positions in the quark_line
	for (uint pos = 0; pos
	< Bas_n_g_minus_2.ca.at(bv).cs.at(ql).ql.size(); pos++)
	// cout << "add_one_gluon( Col_amp, int): Looping over pos " << pos <<endl;

	// Jump over index g_ex by increasing the index g_ex, and all indices above with 1
	if (Bas_n_g_minus_2.ca.at(bv).cs.at(ql).ql.at(pos)
	>= g_ex) {
	// Change the index if >= g_ex
	Cs_new.cs.at(ql).ql.at(pos)++;
	}
	}
	//Col_str Cs_new; //To contain the new basis vector
	//Cs_new=Bas_n_g_minus_2.ca.at(bv);
	Cs_new.cs.push_back(Ql2_ring);
	New_bas = New_bas + Cs_new;
	//cout << "add_one_gluon( Col_amp, int): New state " << Cs_new
	// << endl;
	}

	}
	}
	return New_bas;
	}

	// Function for creating a basis with n_q=n_qbar quarks and n_g gluons.
	// Note: If there are only gluons the anti-cyclic part of the Col_amp is implicit.
	// For example, for (1,2,3) there is also an implicit part -(1,3,2) which is
	// not separately constructed. This is compensated for in scalar_product_matrix and scalar_product_matrix_mem,
	// but not in scalar_product()
	Col_amp Col_functions::create_trace_basis( int n_q, int n_g, int n_loop ) const{
	//cout << "create_trace_basis: Called with n_g=" << n_g << " and n_q=" << n_q << endl;

	// To contain the resulting basis
	Col_amp Basis;

	// First connect all q and qbars in all n_q! ways
	// This gives a basis for the case of n_q quarks and n_qbar anti-quarks
	if( n_q != 0) Basis=connect_quarks( n_q );

	// If there were no quarks
	if( n_q==0 ){
	// There has to be at least two gluons
	if( n_g <= 1 ) {
	std::cerr << "create_trace_basis: For 0 quarks there is no basis with only " << n_g << " gluons" << std::endl;
	throw -1;
	}
	// If 2 or more gluons, build from the 2-gluon basis
	else if( n_g>=2 ){
	Col_str Cs_OnlyState("[(1,2)]");
	Col_amp Ca_tmp;
	Ca_tmp.ca.push_back(Cs_OnlyState);

	Basis= Ca_tmp;
	// For 2 gluons, the work is done
	if( n_g==2 ) return Basis;
	//cout << Basis;
	}
	}

	// Then, add the gluons one at the time

	// If there are only gluons the generation should start from gluon 3,
	// otherwise from 2*n_q+1;
	for( int g_new=std::max(3, 2n_q+1); g_new <= 2n_q+ n_g; g_new++){
	// If only gluons start from the 2-gluon state, so add gluon 3
	//cout << "create_trace_basis: gluon to be added "<< g_new<< endl;
	Basis=add_one_gluon(Basis, n_q ,n_g, g_new, n_loop);
	//add_one_gluon( Col_amp Old_bas, int n_g, int n_q, int g_new)
	}

	// If there are only gluons, to count only physical states each basis state
	// should be combined with its complex conjugate
	// For uniformity, do this for the 2-gluon states as well

	// Normal order the Col_str's
	Basis.normal_order();
	return Basis;

	}

	/*
	// A function for radiating n_g gluons n_run times
	double Col_functions::just_radiate( Col_amp & Ca1, int n_g_rad, double n_run, int n_loop ) const{

	int n_q=Ca1.n_quark();
	int n_g=Ca1.n_gluon();

	// Find the basis
	std::cout << "just_radiate: Constructing the basis for " << n_q
	<< " q qbar pair(s) and " << n_g + n_g_rad << " gluon(s)..."<< std::endl;
	std::cout.flush();
	Col_amp Basis = create_trace_basis(n_q, n_g + n_g_rad, n_loop);

	std::cout << "just_radiate: Constructed the basis: " << Basis << " with " << Basis.ca.size() << " basis tensors \n";
	std::cout.flush();
	// and calculate the scalar product matrix in this basis
	std::cout << "just_radiate: Calculating the scalar product matrix..." << std::endl;
	// Calculating the scalar product matrix and its numerical versions (once and for all)
	// do this both in the leading Nc limit and else
	poly_matr sp_m = scalar_product_matrix(Basis);
	poly_matr sp_m_leading = leading(sp_m);
	dmatr sp_m_num = double_num(sp_m);
	dmatr sp_m_num_leading = double_num(sp_m_leading);
	std::cout << sp_m;
	//cout << sp_m_leading;

	// For containing the scalar product
	cnum sp_num, sp_num_leading;

	// For containing the average ratio of amp^2 for Nc=3 and Nc=infinity
	double ratio_average=0;

	// Col_amp to start the radiation from
	Col_amp Ca_initial = Ca1;

	// The initial number of partons
	int n_initial = 2 * n_q + n_g;

	std::cout << "just_radiate: Running the shower... " << std::endl;

	// Make n different showers
	for (int n = 1; n <= n_run; n++) {
	std::cout << "just_radiate: Calculating amplitude... " << std::endl;

	// Reset n_parton and the starting Col_amp
	int n_parton = n_initial;
	Ca1 = Ca_initial;

	// Radiate n_g partons
	// Radiate parton p, where p=n_parton+1 ... n_initial+n_g_rad
	for (int p = n_parton + 1; p <= n_initial + n_g_rad; p++) {

	Col_amp Ca_new;
	// Add contributions where p is radiated from all partons (one at the time)
	for (int i = 1; i <= n_parton; i++) {
	// Radiate parton p from parton i Ca
	Col_amp part_i = emit_gluon(Ca1, i, p);


	// Multiply with a random number to mimic amplitude from other DOF
	double r=1.0;
	r=rand()/float(RAND_MAX)-0.5;
	part_i=part_i*r;

	// Append to new Ca (the Ca after emission of one more parton)
	Ca_new.append(part_i.ca);
	}
	// Now the number of partons have increased with 1
	n_parton++;

	// Simplify the resulting amplitude (turning off makes things slow, but tests scalar prod)
	Ca_new.simplify();
	Ca1 = Ca_new;
	//cout << "just_radiate: Ca1" << Ca1 << endl;
	}

	// cout << "just_radiate: Ca1" << Ca1 << endl;
	// Calculate the scalar product
	std::cout << "just_radiate: Calculating \|A\|^2..." << std::endl;
	//sp = scalar_product(Ca1, Ca1, Basis, sp_m);

	sp_num = scalar_product_cnum_num(Ca1, Ca1, Basis, sp_m_num);
	sp_num_leading = scalar_product_cnum_num(Ca1, Ca1, Basis, sp_m_num_leading);

	// Write out result the first few times
	if (n <= 2) {
	std::cout << "Resulting Col_amp" << Ca1 << std::endl;
	std::cout << "with scalar prod " << sp_num << std::endl;
	std::cout << "In the infinity limit " << sp_num_leading << std::endl;
	}

	std::cout << "Run " << n << " The ratio Nc=Infinity / Nc=3 " <<
	sp_num_leading/sp_num << std::endl;
	ratio_average = ratio_average + real(sp_num_leading / sp_num ) / n_run;

	}

	return ratio_average;
	}
	*/



	// Function for creating the name of the file containing the scalar product matrix
	// The argument leading should be true if it is a leading matrix that is to be saved.
	// Function for creating the name of the file containing the radiation amplitude matrix.
	// The argument leading should be true if it is a leading matrix that is to be saved.
	std::string Col_functions::ram_file_name( const Col_amp & Basis, bool leading ) const{

	// First construct filename
	std::ostringstream ss;
	std::string filename;
	ss << "ColorResults";
	ss << '/';
	ss << "RAM_q";
	int nq=Basis.at(0).n_quark();
	ss << nq;
	ss << "_g";
	int ng=Basis.at(0).n_gluon();
	ss << ng;
	if ( leading ) ss << "_l";
	if ( leading_CF ) ss << "_CFl";
	ss << "_Nc";
	ss << Nc;
	ss << ".dat";
	filename=ss.str();
	return filename;
	}


	void Col_functions::write_out_spm( const dmatr & matr, const char* filename ) const{
	std::ofstream outfile(filename);
	outfile << matr;
	}

	void Col_functions::write_out_spm( const Poly_matr & matr, const Col_amp & Basis, bool leading ) const{

	// Creating file name (part)
	std::string filename = spm_file_name(Basis, leading);

	// Then, write to the file
	std::ofstream outfile(filename.c_str());

	outfile << matr;
	}

	// Function for writing out a numerical scalar product matrix to the default filename and file path
	// The filename is constructed using info of the basis
	cnum Col_functions::scalar_product_num_diagonal( const Col_amp & Ca1, const Col_amp & Ca2, const dmatr & sp_matr ) const{

	// To contain the result
	cnum res=0;
	uint basis_size=Ca1.ca.size(); // Both Ca's and sp_matr should have same length
	time_t t1, t2, t3;
	t1=clock( );

	// Make ordinary vectors of numerical information in Polynomials, to speed up multiplication
	cvec v1, v2;
	v1.reserve(basis_size);
	for (uint i=0; i< basis_size; i++){
	// added conjugate of first arg 130524
	v1.push_back( conj( cnum_num( Ca1.ca.at(i).Poly ) ) );
	};
	if( &Ca1==&Ca2 ) v2=v1;
	else{
	v2.reserve(basis_size);
	for (uint i=0; i< basis_size; i++) v2.push_back( cnum_num( Ca1.ca.at(i).Poly ) );
	}

	t2=clock();
	std::cout << "scalar_product_num_diagonal(Ca1,Ca2,dmatr): Decomposition done at time: " << t2-t1 << std::endl;

	// Then add contributions
	for (uint m1=0; m1< basis_size; m1++){
	res=res+ v1.at(m1) v2.at(m1) sp_matr.at(m1).at(m1);
	}
	t3=clock();
	std::cout << "scalar_product_num_diagonal(Ca1,Ca2,dmatr): Matrix multiplication done at time: " << t3-t2 << std::endl;
	return res;

	}

	void Col_functions::write_out_spm( const dmatr & matr, const Col_amp & Basis, bool leading ) const{
	//cout << "write_out_spm(dmatr, Ca, bool): Incoming matrix \n" << matr;

	// Creating file name (part)
	std::string filename = spm_file_name(Basis, leading);

	// Then, write to the file
	std::ofstream outfile(filename.c_str());

	outfile << matr;
	}



	// Function for calculating the weight < M \| T^(i) T^(j) \| M > / (< M \| M > (c_f or c_a) ),
	// used in the Sudakov for generating the next k_T.
	// The Ca should be \| M >, and i and j are the partons evolved in the emission.
	// The factor (c_f or c_a) is present to compensate for factor built in in kinematics part
	cnum Col_functions::color_shower_weight( Col_amp Ca, int i, int j, int g_new ) const{

	// The amplitudes after emission
	Col_amp Cai = emit_gluon(Ca ,i, g_new);
	Col_amp Caj = emit_gluon(Ca ,j, g_new);

	cnum numerator=0;
	cnum denominator=1.0;

	int n_q=Ca.n_quark();
	int n_g=Ca.n_gluon();

	// If gluons only, take care of implicit scalar product
	if( Ca.gluons_only() ){
	// The numerator
	int sign = ((n_g+1) % 2 ? -1 : 1);
	Col_amp Caj_conj=Caj;
	Caj_conj.conjugate();
	numerator = cnum_num( 2scalar_product (Cai, Caj) +2sign*scalar_product(Cai,Caj_conj) );

	// Before the emission there is one parton less, so the sign is opposite
	sign = (n_g % 2 ? -1 : 1);
	Col_amp Ca_conj=Ca;
	Ca_conj.conjugate();
	denominator = cnum_num( 2scalar_product(Ca, Ca) +2sign*scalar_product (Ca, Ca_conj ) );

	} else if(n_q==1){
	numerator = cnum_num( scalar_product (Cai, Caj) );
	denominator = cnum_num( scalar_product (Ca, Ca) );
	}

	// The CF or CA factor
	if ( Ca.ca.at(0).find_kind(i) == "g" ) denominator=denominator*Nc;
	else denominator=denominator*CF;

	return numerator/denominator;
	}

	// Function for calculating the radiation amplitudes associated with emissions from
	// the various partons, starting with n_q quarks and n_g gluons.
	// The content is C^m T^(i) T^(j) C^n as a matrix of matrices.
	// The factors associated with i and j stand in the same sub-matrix.
	// For example, the first n_basis*n_basis block corresponds to the coherent
	// emission from parton 1 and parton 1. The block to the right of this
	// corresponds to coherent emission of parton 1 and parton 2 etc.
	// Blocks further down in the matrix corresponds to first parton > 1.
	// TODO test leading terms make work for nq!=1
	dmatr Col_functions::radiation_amplitue_matrix( uint n_q, uint n_g, int n_loop ) const{
	std::cout << "radiation_amplitue_matrix: calculating radiation amplitue matrix" << std::endl;
	// Total number of partons, for each q there is a qbar
	uint n_partons=2*n_q+n_g;
	// The gluon to be emitted
	int g_new=n_partons+1;

	// Make sure only used at tree level
	if (! (n_loop==0 && (n_q==1 or n_q==0)) ) {
	std::cerr << "radiation_amplitue_matrix: Function so far only intended for tree level processes with one open Quark_line." << std::endl;
	throw -1;}

	// The basis before emission
	std::cout << "radiation_amplitue_matrix: creating basis for " << n_q << " quark(s) (and as many anti-quarks) and "
	<< n_g << " gluons..." << std::endl;
	Col_amp Basis=create_trace_basis(n_q, n_g, n_loop);
	uint basis_size=Basis.size();
	std::cout << "radiation_amplitue_matrix: Constructed basis with " << Basis.size() << " basis vectors." << std::endl;

	// For remembering already calculated topologies
	std::map< std::string, Polynomial> mem_map;

	// To contain the resulting color weight matrices _l is the leading version
	dmatr matr;
	dmatr matr_l;
	matr.reserve( basis_size*n_partons );
	matr_l.reserve( basis_size*n_partons );

	// Initialize matrices to 0
	for (uint i = 0; i < n_partons*basis_size; i++) {
	std::vector<double> rowi, rowi_l;// Create the row i
	for (uint j = 0; j < n_partons*basis_size; j++) {
	rowi.push_back(0); // Create entry j=0
	rowi_l.push_back(0); // Create entry j=0
	}
	matr.push_back(rowi);
	matr_l.push_back(rowi);
	}
	std::cout << "radiation_amplitue_matrix: Calculating 'topologies'..." << std::endl;
	// If there is only one q qbar pair get all topologies by allowing the q and qbar (1 and 2)
	// and 2 different gluons to radiate
	int partons_needed=n_partons;
	if(n_q==1) partons_needed=4;
	if(n_q==0) partons_needed=4; // Gluon 4 can stand in any position
	// Loop over possible emitters i and j
	for (int i = 1; i <= partons_needed; i++) {
	std::cout << "radiation_amplitue_matrix: Calculating 'topologies' for emitter " << i << std::endl;

	for (int j = 1; j <= partons_needed; j++) {
	// Loop over basis vectors (in the relevant basis before emission)
	for (uint m = 0; m < basis_size; m++) {
	for (uint n = 0; n < basis_size; n++) {
	//cout << "i " << i << " j " << j << " m " << m << " n " << n << endl;

	// Calculate element
	// The amplitudes after emission
	Col_amp Cai = emit_gluon(Basis.at(m), i, g_new);
	Col_amp Caj = emit_gluon(Basis.at(n), j, g_new);

	// Rename indices, and make string of new Col_amps, to use in map
	std::pair<Col_amp, Col_amp> Css_new = rename_indices(Cai, Caj);
	Col_amp Ca1 = Css_new.first;
	Col_amp Ca2 = Css_new.second;
	std::ostringstream Ca_string;
	Ca_string << Ca1 << Ca2;
	//cout << double_num(scalar_product(Caj, Cai)) << " ";

	// Add result to map
	if (!(mem_map.count(Ca_string.str()) > 0)) {
	//cout << "adding term for" << endl << Ca_string.str();
	Polynomial p; // to contain result
	if(n_q==1) p = scalar_product (Cai, Caj);
	else if(n_q==0) {
	// Add implicit gluon state
	int sign = ((n_g+1) % 2 ? -1 : 1);
	Col_amp Cai_conj=Cai;
	Cai_conj.conjugate();
	p = 2scalar_product (Cai, Caj) + 2sign*scalar_product( Cai_conj, Caj ) ;
	}

	p.simplify();
	mem_map.insert(make_pair(Ca_string.str(), p));
	}

	} // End of m-loop
	} // End of n-loop
	} // End of j-loop
	} // End of i-loop

	std::cout << "radiation_amplitue_matrix: Calculating numerical version..." << std::endl;
	std::map< std::string, double> mem_map_num=double_num(mem_map);

	// Taking the leading part for each "topology" and making numerical maps
	std::cout << "radiation_amplitue_matrix: Calculating leading numerical version..." << std::endl;
	std::map< std::string, double> mem_map_l_num=double_num(leading(mem_map));


	// After this the mem_map is not needed anymore
	mem_map.clear();
	std::cout << "radiation_amplitue_matrix: Calculating matrix entries..." << std::endl;

	// Loop over possible emitters i and j
	for (uint i = 1; i <= n_partons; i++) {
	std::cout << "radiation_amplitue_matrix: Calculating matrix entries for emitter " << i << std::endl;
	for (uint j = 1; j <= n_partons; j++) {
	// Loop over basis vectors (in the relevant basis before emission)
	for (uint m = 0; m < basis_size; m++) {
	for (uint n = 0; n < basis_size; n++) {

	// The amplitudes after emission
	Col_amp Cai = emit_gluon(Basis.at(m), i, g_new);
	Col_amp Caj = emit_gluon(Basis.at(n), j, g_new);

	// Rename indices, and make string of new Col_amps, to use in map
	std::pair<Col_amp, Col_amp> Css_new = rename_indices(Cai, Caj);
	Col_amp Ca1 = Css_new.first;
	Col_amp Ca2 = Css_new.second;
	std::ostringstream Ca_string;
	Ca_string << Ca1 << Ca2;

	// Find the value using the map
	double entry = mem_map_num[Ca_string.str()];
	double entry_l = mem_map_l_num[Ca_string.str()];

	// Insert element in right place
	uint row=basis_size*(i-1)+m;
	uint column= basis_size*(j-1)+n;

	matr.at(row).at(column)=entry;
	matr_l.at(row).at(column)=entry_l;
	}
	}
	}
	}

	// Verifying that the matrix is symmetric to accuracy
	std::cout << "radiation_amplitue_matrix: Checking symmetry... ";
	bool sym = check_symmetry(matr);
	if (!sym)
	throw -1;

	// Also save the numerical versions to files for future usage
	std::cout << "radiation_amplitue_matrix: Writing out color radiation amplitude matrix... " << std::endl;
	std::ofstream outfile(ram_file_name(Basis, false).c_str());
	outfile << matr;
	//matr.clear();
	std::cout << "done." << std::endl;
	std::cout << "radiation_amplitue_matrix: Writing out leading contributions to radiation amplitude matrix... " << std::endl;
	std::ofstream outfile_l(ram_file_name(Basis, true).c_str());
	outfile_l << matr_l;
	std::cout << "done." << std::endl;
	matr_l.clear();
	std::cout << "radiation_amplitue_matrix: Returning matrix... " << std::endl;

	return matr;
	}


	// Function for calculating the weight < M \| T^(i) T^(j) \| M > / (< M \| M > (c_f or c_a) ),
	// used in the Sudakov, for generating the next k_T.
	// The amp should be \| M >, and i and j are the partons involved in the emission.
	// The index g_new denotes the parton to be emitted and spm the scalar product matrix used
	// for numerical evaluation (before emission).
	// The factor (c_f or c_a) is present to compensate for factor built in in kinematics part.
	// In the momentum code i is the emitter and j is the spectator.
	// The matrix radamps contains the amplitudes associated with emission,
	// and can be calculated using radiation_amplitue_matrix.
	// The number of existing quarks is n_q and the number of gluons is n_g.
	cnum Col_functions::color_shower_weight( const cvec & amp, int i, int j, int g_new, const dmatr & radamps, const dmatr & spm_matr, int n_q, int n_g, int n_loop ) const{


	// Make sure only used at tree level
	if (! (n_loop==0 && (n_q==1 or n_q==0)) ) {
	std::cerr << "color_shower_weight: Function so far only intended for tree level processes with one quark line." << std::endl;
	throw -1;}

	//int n_partons=2*n_q+n_g;

	// Calculate the size of the basis before emission
	//uint basis_size=radamps.size()/n_partons;
	uint basis_size=spm_matr.size();


	// To contain the numerator and denominator
	cnum numerator=0;
	cnum denominator=1.0;

	// The denominator <M\| SPM \|M>
	// This should work also in the gluons only case, as states implicit in spm_matr and in amp
	denominator = scalar_product_num_fast(amp, amp, spm_matr);

	// The CF or CA factor (1 is a q 2 is a qbar)
	if ( (i==1 or i==2) && n_q==1) denominator=denominator*CF;
	else denominator=denominator*Nc;


	// To find the numerator loop over color structures
	// Loop over basis vectors in \| M >
	for (uint m = 0; m < basis_size; m++) {

	// Loop over basis vectors in < M \|
	for (uint n = 0; n < basis_size; n++) {

	uint row=basis_size*(i-1)+m;
	uint column= basis_size*(j-1)+n;
	// This should work also in the gluons only case, as states implicit in rad_amps and in amp
	numerator += amp.at(m)radamps.at(row).at(column)amp.at(n);
	}
	}

	return numerator/denominator;
	}

	// The factorial of an int
	int Col_functions::factorial (int i) const{
	if(i<0) {
	std::cerr << "Col_functions: factorial: intended for int >=0, argument was " << i << std::endl;
	std::cerr.flush();
	throw -1;
	}
	if (i==0) return 1;
	return factorial(i-1)*i; // Recursive call
	}

	void Col_functions::write_out_dvec( const dvec & dv, std::string filename ) const {

	std::ofstream outfile(filename.c_str() );
	outfile << dv;
	}

	dmatr Col_functions::read_in_dmatr( std::string filename ) const {


	// Read in file
	std::ifstream fin(filename.c_str());

	// Check that file exists
	if( !fin ){
	std::cerr << "Col_functions::read_in_dmatr: The file "
	<< filename << " could not be opened." << std::endl;
	throw filename;

	}

	// Copy info from file to string
	std::string str((std::istreambuf_iterator<char>(fin)), std::istreambuf_iterator<char>());

	// Skip lines starting with #
	while(str.at(0)== '#'){
	while (str.at(0) != '\n'){
	str.erase(str.begin());
	}
	// erase endl sign(s)
	while(str.at(0)== '\n'){
	str.erase(str.begin());
	}
	//std::cout << "Col_basis::read_in_dmatr: erased comment, now sting is " << str << std::endl;
	}

	// First char in file should be '{'
	if (str.at(0) != '{') {
	std::cerr
	<< "Col_functions::read_in_dmatr: First char in matrix data file after comments should be '{', it was: "
	<< str.at(0) << std::endl;
	throw -1;
	}


	// Check that only allowed signs
	uint j = 0;
	while (j < str.size()) {

	if (!(str.at(j) == '+' or str.at(j) == '-' or str.at(j) == '.'
	or str.at(j) == '{' or str.at(j) == '}' or str.at(j) == '\n'
	or str.at(j) == ',' or str.at(j) == ' ' or str.at(j) == '0'
	or str.at(j) == '1' or str.at(j) == '2' or str.at(j) == '3'
	or str.at(j) == '4' or str.at(j) == '5' or str.at(j) == '6'
	or str.at(j) == '7' or str.at(j) == '8' or str.at(j) == '9')) {
	std::cerr
	<< "Col_functions::read_in_dmatr: A disallowed sign encountered in string for dmatr: "
	<< str.at(j) << ", in file " << filename << std::endl;
	std::cerr << "Col_functions::read_in_dmatr expects a numerical matrix." << std::endl;
	throw str;
	}
	j++;
	};


	// Row to contain numbers
	dvec row;

	// To contain matrix of scalar products
	dmatr matr;

	// Read the string, starting from 0th element
	unsigned int i = 0;
	while (i < str.size() - 2) {
	i += 1;

	// We may have to skip some chars
	while (i< str.size()-2 &&(str.at(i) == ',' or str.at(i) == '}' or str.at(i) == ' ' or str.at(i) == '\n' or str.at(i) == ' ' or str.at(i) == '{') ) i++;

	// String to make a number of, and double to contain number
	std::string num_str;
	num_str.clear();
	double num;

	// Keep reading the number while not ',' or '}'
	while ( i< str.size()-2 && (str.at(i) != ',' && str.at(i) != '}') )
	{
	num_str.push_back(str.at(i));
	i++;
	//std::cout << "Col_basis::read_in_dmatr at " << i << " "
	// << num_str << std::endl;

	}
	// now, at(i), there is either , or }
	//std::cout << str.at(i);

	// num_str contains the string to make a number of
	std::istringstream parton_str_st( num_str );
	parton_str_st >> num;
	//std::cout << "Col_basis::read_in_dmatr num at " << i << " "
	// << num << std::endl;


	// Add number to vector
	row.push_back(num);

	// If we have a new row
	if( i< str.size()-2 && str.at(i)=='}'){
	// Save row in matrix, and empty row
	matr.push_back(row);
	row.clear();

	//std::cout << "Intermediate matr : " << matr << std::endl;
	//std::cout << str.at(i);

	// We may have to skip some chars
	while (i< str.size()-2 &&(str.at(i) == ',' or str.at(i) == '}' or str.at(i) == ' ' or str.at(i) == '\n' ) ) {
	i++;
	}
	}

	// Otherwise just keep on reading the next number in row
	};

	return matr;
	}

	dvec Col_functions::read_in_dvec( std::string filename ) const {


	// Read in file
	std::ifstream fin(filename.c_str());

	// Check that file exists
	if( !fin ){
	std::cerr << "Col_functions::read_in_dvec: The file "
	<< filename << " could not be opened." << std::endl;
	throw filename;

	}

	// Copy info from file to string
	std::string str((std::istreambuf_iterator<char>(fin)), std::istreambuf_iterator<char>());

	// Skip lines starting with #
	while(str.at(0)== '#'){
	while (str.at(0) != '\n'){
	str.erase(str.begin());
	}
	// erase endl sign(s)
	while(str.at(0)== '\n'){
	str.erase(str.begin());
	}
	//std::cout << "Col_basis::read_in_dmatr: erased comment, now sting is " << str << std::endl;
	}

	// First char in file should be '{'
	if (str.at(0) != '{') {
	std::cerr
	<< "Col_functions::read_in_dvec: First char in matrix data file after comments should be '{', it was: "
	<< str.at(0) << std::endl;
	throw -1;
	}


	// Check that only allowed signs
	uint j = 0;
	while (j < str.size()) {

	if (!(str.at(j) == '+' or str.at(j) == '-' or str.at(j) == '.'
	or str.at(j) == '{' or str.at(j) == '}' or str.at(j) == '\n'
	or str.at(j) == ',' or str.at(j) == ' ' or str.at(j) == '0'
	or str.at(j) == '1' or str.at(j) == '2' or str.at(j) == '3'
	or str.at(j) == '4' or str.at(j) == '5' or str.at(j) == '6'
	or str.at(j) == '7' or str.at(j) == '8' or str.at(j) == '9')) {
	std::cerr
	<< "Col_functions::read_in_dvec: A disallowed sign encountered in string for dmatr: "
	<< str.at(j) << ", in file " << filename << std::endl;
	std::cerr << "Col_functions::read_in_dvec expects a numerical matrix." << std::endl;
	throw str;
	}
	j++;
	};


	// Row to contain numbers
	dvec row;

	// To contain matrix of scalar products
	dmatr matr;

	// Read the string, starting from 0th element
	unsigned int i = 0;
	while (i < str.size() - 1) {
	i += 1;

	// We may have to skip some chars
	while (i< str.size()-2 &&(str.at(i) == ',' or str.at(i) == '}' or str.at(i) == ' ' or str.at(i) == '\n' or str.at(i) == ' ' or str.at(i) == '{') ) i++;

	// String to make a number of, and double to contain number
	std::string num_str;
	num_str.clear();
	double num;

	// Keep reading the number while not ',' or '}'
	while ( i< str.size()-1 && (str.at(i) != ',' && str.at(i) != '}') )
	{
	num_str.push_back(str.at(i));
	i++;
	//std::cout << "Col_basis::read_in_dmatr at " << i << " "
	// << num_str << std::endl;

	}
	// now, at(i), there is either , or }
	//std::cout << str.at(i);

	// num_str contains the string to make a number of
	std::istringstream num_str_st( num_str );
	num_str_st >> num;
	//std::cout << "Col_basis::read_in_dmatr num at " << i << " "
	// << num << std::endl;


	// Add number to vector
	row.push_back(num);

	// Skip signs in end and make suere not to enter loop one extra time
	while( i<str.size() and ( str.at(i)==' ' or str.at(i)=='\n' or str.at(i)=='}' ) ) i++;

	};

	return row;

	}

	void Col_functions::write_out_dmatr( const dmatr & matr, std::string filename ) const {
	std::ofstream outfile(filename.c_str());
	outfile << matr;
	}

	// Function for finding out the number of the vector with only one quark_line
	// with the color information of ql.
	// The quark_line should be normal ordered before usage.
	int Col_functions::quark_line_number( int n_q, int n_g, const quark_line & ql ) const{

	if( !( ql.at(0)==1 && n_q==1 && ql.at( ql.size()-1 )==2 ) && !( ql.at(0)==1 && n_q==0) ){
	std::cerr << "quark_line_number: Expects a normal ordered quark_line with 1 as the first entry "
	<< " (and 2 the last entry if only gluons)."
	<< std::endl << "Received " << ql << std::endl;
	throw -1;
	}
	if( n_q==0 ){
	std::cerr << "quark_line_number: Does not yet work for gluons only." << std::endl;
	throw -1;
	}

	// To contain the number of the vector
	int tens_num=0;

	// If there is only one vector
	if( (n_q==0 && n_g==2) or (n_q==1 && n_g==1) )
	return 0;

	// Loop over positions and count up tenors number
	// pos 0 should contain 1 and is irrelevant, pos last is similarly irrelevant if quark_line has quarks.
	for(int pos=1; pos <= static_cast<int>(ql.size())-1-n_q; pos++){
	// Look among already found and see how many where smaller
	// The effective parton number shoudl be counted down with this number
	// Alternatively one can look to the right and see how many where larger
	// Thus what counts is the number of smaller integers to the right
	int parton=ql.at(pos);
	for(int p=1; p< pos; p++)
	if( ql.at(p) < ql.at(pos)) parton=parton-1;
	//cout << "factorial " << factorial(n_g+n_q-1-pos) << endl;
	//cout << "other part " << (parton-1-n_q-1) << endl;
	tens_num=tens_num+factorial(n_g+n_q-1-pos)*(parton-1-n_q-1);
	}

	return (tens_num);
	}


	// Function for locating the parton parton in the normal ordered basis,
	// given the number of the vector vecnum_num, and the number of quarks and gluons in the basis
	std::pair<int,int> Col_functions::find_parton( int parton, int vecnum_num, int n_q, int n_g, int n_loop ) const{
	if( !(n_q==1 or n_q==0 or n_q==2) or n_loop >0 ){
	std::cerr << "find_parton: Function only intended for special case of 0-2 q qbar pairs, and a tree level basis." << std::endl;
	throw -1;
	}
	if( parton==1 && n_q == 1) return std::make_pair( 0, 0 ); // q=1 or g=1 has position 0 if one ql
	if( parton==1 && n_q == 0) return std::make_pair( 0, 0 ); // q=1 or g=1 has position 0 if one ql
	if( parton==2 && n_q==1 ) return std::make_pair( 0, n_g+1 ); // qbar has last pos

	// If n_q=2 the color structure is of type {q1, g1....gm,q2}{q3,gm+1...gn,q4}
	// This is almost as in the case of only one Ql {q1, g1....gm,gm+1...gn,q4}
	// with a breakpoint "q2}{q3" inserted or q4}{q3.
	// To find the tensor number, first locate the breakpoint.
	// This can be done by noting that there are 2*n_g! options for each breakpoints,
	// where the factor 2 comes from the 2 options for the first anti-quark
	// we thus have
	int break_after=n_g;
	if (n_q == 2) {
	// There are n_q!n_q!n_g! assignments for each break
	break_after = break_after-(vecnum_num) / (22factorial(n_g));
	vecnum_num= vecnum_num -(n_g-break_after)(22*factorial(n_g));
	//cout << "find_parton: break_after " << break_after << endl;
	//cout << "find_parton: vecnum_num after break compensation " << vecnum_num << endl;

	// Check if first parton is 1 or 3, there are 2 options for each from the qbar being 2 or 4
	bool first3=false;
	if ( ( vecnum_num / ( 2*factorial(n_g)) ) ) first3=true;
	//cout << "The first is 3? "<< first3;

	vecnum_num= vecnum_num % (2*factorial(n_g));

	// Special cases of asked for quark
	if(parton == 1){
	if(!first3) return std::make_pair(0,0);
	else return std::make_pair(1,0);
	}
	if(parton == 3){
	if(first3) return std::make_pair(0,0);
	else return std::make_pair(1,0);
	}
	//cout << "find_parton: vecnum_num after first quark compensation " << vecnum_num << endl;
	}


	// The number of partons gluons with lower parton number than the parton
	// (and possibly standing to the right, hence -1)
	int lower=parton-2-n_q;
	//if (n_q==2 ) lower=parton -2*n_q-1;
	if ( n_q==2 ) lower=parton -2*n_q-1; // quarks don't count


	// To contain the number of gluons standing to the right of the parton under consideration
	// and being smaller than 2, noting to the right is smaller than 2, hence 0
	// needed for special case of gluons only
	int lower2=0;

	// Number of partons to the right, when checking the various positions
	int n_r=0;
	// For quark lines(s) all gluons except the one under consideration matter
	if ( n_q>0 ) n_r=n_g-1;
	// For a gluon line we have to subtract 1 extra as we start read at place 1 anyway
	else if ( n_q==0 ) n_r=n_g-2;

	// Was the gluon 2 found or not, needed for special case of gluons only
	bool found_2=false;

	// The resulting quark_line number (can change only in case 2 q qbar pairs)
	int ql_num=0;

	// Loop over possible places of the gluon, start with checking place 1
	// n_l_r contains the number of partons lower than the parton at the position under consideration
	int n_l_r=0; // To make compiler not warn
	int looped_over=0;
	while(true){
	if( n_q==1 ) {
	n_l_r=vecnum_num/factorial(n_r); // The result of integer division
	vecnum_num=vecnum_num%factorial(n_r); // The rest after integer division
	}

	if( n_q==2 ) {

	// If both qbars to right extra factor 2!
	if( break_after > looped_over) {
	//cout << "break_after > looped_over=" << looped_over << endl;

	n_l_r=vecnum_num/(factorial(n_r)*2); // The result of integer division
	vecnum_num=vecnum_num%(factorial(n_r)*2); // The rest after integer division
	}
	else{
	// Check if it is time to compensate for first qbar and second q
	if ( break_after==looped_over ){
	//cout << "break_after == looped_over=" << looped_over << endl;
	//cout << "qbar compensation starting with vector number " << vecnum_num << endl;
	//cout.flush();

	// If we have already looped over all gluons
	if(n_r==-1) n_r=0;
	//cout << "find_parton: n_r " << n_r << endl;
	//cout.flush();
	// Special cases if asked for first qbar = 2 or 4
	if( parton==2){
	// All gluons which are not looped over are to right
	if (vecnum_num/(factorial(n_r+1))==0) return std::make_pair(0,break_after+1);
	else return std::make_pair(1,n_g-break_after+1);
	}
	if( parton==4){
	if (vecnum_num/(factorial(n_r+1))==1) return std::make_pair(0,break_after+1);
	else return std::make_pair(1,n_g-break_after+1);
	}
	//cout << "find_parton: vecnum_num before qbar compensation " << vecnum_num << endl;
	//cout << "find_parton: n_r " << n_r << endl;


	vecnum_num=vecnum_num%(factorial(n_r+1)); // compensated for first qbar


	//cout << "find_parton: vecnum_num after first qbar compensation " << vecnum_num << endl;
	// After looking at first qbar the ql_num should be 1
	ql_num=1;
	}
	//else cout << "break_after < looped_over=" << looped_over << endl;



	// Otherwise the situation is as for gluons only
	n_l_r=vecnum_num/(factorial(n_r)); // The result of integer division
	vecnum_num=vecnum_num%factorial(n_r); // The rest after integer division
	}
	//cout << "n_l_r " << n_l_r << endl;
	}
	else if (n_q==0){
	//cout << "found 2? " << found_2 << endl;
	if( found_2 ){
	//cout << "New n_r " << n_r << endl;
	n_l_r=vecnum_num/factorial(n_r);
	vecnum_num=vecnum_num%factorial(n_r);
	//cout << "n_l_r " << n_l_r << endl;
	}
	else{
	// First check if 2 was at the position
	//cout << "Before finding 2 n_r " << n_r << endl;
	//cout << "Before finding 2 vecnum_num " << vecnum_num << endl;
	n_l_r=vecnum_num/(factorial(n_r)); // If 2 was at the place, then all orders matter

	//cout << "Before finding 2 n_l_r " << n_l_r << endl;

	// If 2 was found
	if( n_l_r==lower2 ) {
	found_2=true;
	vecnum_num=vecnum_num%(factorial(n_r)); }
	else{ // Recalculate knowing that we did not find 2
	n_l_r=vecnum_num/(factorial(n_r)/2); // There is only one relative order of gluon 2 and 3
	vecnum_num=vecnum_num%(factorial(n_r)/2); // The rest after integer division
	}

	}
	}
	if( n_l_r==lower ) break; // Number smaller to right = number smaller than parton to right
	if( n_l_r < lower ) lower--;
	n_r--;
	looped_over++; // needed to compare with break in 2q case
	//cout << "n_l_r " << n_l_r << endl;
	}
	//if( n_q==1) return make_pair( 0, n_g-n_r );
	//else if( n_q==0 ) return make_pair( 0, n_g-n_r-1 );

	// If we have 2 quark lines we have to first figure out the quark line number
	//int ql_num=n_g-break_after;

	// and then the number in the quark_line
	int pos_num=n_g-n_r; // n_q==1 case
	if( n_q==0) pos_num=n_g-n_r-1;
	else if( n_q==2 && ql_num==0) pos_num=n_g-n_r;
	else if(n_q==2 && ql_num==1) pos_num=n_g-break_after-n_r;

	return std::make_pair( ql_num, pos_num );
	}

	// Function for finding the new vector number in the basis with n_p+1 partons
	// after inserting a new parton with larger parton number at the place place.
	// The old vector has number oldouble_num, and there were, before emission
	// n_q quarks (+ n_q anti-quarks) and n_g gluons, i.e. n_p=2 n_q+ n_g.
	// Function so far only intended for tree-level processes with at most 2 qqbar pairs.
	// Function explicitly tested for initial states with 2qqbar pairs and up to 5 gluons,
	// 1qqbar pair and up to 7 gluons, 0 qqbar pairs and up to 8 gluons.
	int Col_functions::new_vector_number( int oldouble_num, std::pair<int,int> place, int n_q, int n_g, int n_loop ) const{
	//cout << "new_vector_number: inserting at place (" << place.first <<", " << place.second << ") for n_g=" << n_g << " and n_q=" << n_q<< endl;
	std::cout.flush();
	if( !( n_q==0 or n_q==1 or n_q==2 ) ){
	std::cerr << "new_vector_number: Function so far only intended for special case of 1 quark line or 2 open quark lines." << std::endl;
	std::cerr.flush();
	throw -1;
	}
	if( n_q==2 && n_g>5 ){
	std::cerr << "new_vector_number: WARNING: function only tested for 2 q qbar pairs and at most 5 gluons." << std::endl;
	}
	if( n_q==1 && n_g>7 ){
	std::cerr << "new_vector_number: WARNING: function only tested for 1 q qbar pair and at most 7 gluons." << std::endl;
	}
	if( n_q==0 && (n_g>8 or n_g==1) ){
	std::cerr << "new_vector_number: WARNING: function only tested 2-8 gluons when no quarks." << std::endl;
	}

	if(place.second == 0 && n_q>0 ){
	std::cerr << "qber: Cannot insert gluon at place 0, reserved for quark." << std::endl;
	std::cerr.flush();
	throw -1;
	} else if(place.second==n_g+2 && n_q==1){
	std::cerr << "new_vector_number: Cannot insert gluon at place n_g+2, reserved for qbar." << std::endl;
	std::cerr.flush();
	throw -1;}
	else if(n_q==2 && (0 > place.second or place.second > n_g+1+1) ){
	std::cerr << "new_vector_number: Cannot insert gluon outside range 1...n_g ." << std::endl;
	std::cerr.flush();
	throw -1;
	}
	else if( 0 > place.second or place.second > n_g+1+n_q){
	std::cerr << "new_vector_number: Cannot insert gluon outside range 1...n_g ." << std::endl;
	std::cerr.flush();
	throw -1;
	}
	// If asked to insert a gluon at place 0, it should be inserted beyond the last place instead
	if( n_q==0 && place.second==0 ){
	place.second=n_g;
	//cerr << "new_vector_number: Cannot insert gluon at place 0 inserting at last place " << place.second << " instead." << endl;
	//cerr.flush();

	}

	// In the gluons only case we need to keep track of the next position relative to gluon 2
	int place2=0;

	int fnext=oldouble_num; // To contain the remaining part of tensor number
	// If there are 2 quarks, there is an effective scalar place in {1,g1...gm,2}{3,gm+1...4}
	// obtained by crossing 2}{3 out. This place is given by
	// (place in the first ql) if it is in the first ql
	// (number of gluons in first ql) + place in second if it's in the second ql

	// If n_q=2 the color structure is of type {q1, g1....gm,q2}{q3,gm+1...gn,q4}
	// This is almost as in the case of only one Ql {q1, g1....gm,gm+1...gn,q4}
	// with a breakpoint "q2}{q3" inserted or q4}{q3.
	// To find the tensor number, first locate the breakpoint.
	// This can be done by noting that there are 22n_g! options for each breakpoints,
	// where the factor 2*2 comes from the options for the quarks
	// we thus have
	int break_after=0;
	int fsplit=0; // Contribution from where the ql is split
	int fq1=0; // Contribution from first quark
	if (n_q == 2) {

	// To contain the number of gluon in the first ql
	// First the number of gluons in 2nd ql
	// For each gluon in 2nd ql there was
	// 2 (from n_q)2 (from n_qbar)factorial(n_g) (from the gluons)
	// tensors with lower numbers from other breaks
	// Note that for all other splitting the factor from the first q is relevant
	break_after = oldouble_num /( 22factorial(n_g) ) ;


	fnext= oldouble_num % (22factorial(n_g));
	// The factor from the split
	fsplit=oldouble_num-fnext;
	break_after=n_g-break_after;
	//cout << "new_vector_number: break_after " << break_after << endl;
	//cout << "new_vector_number: fsplit " << fsplit << endl;
	//cout << "new_vector_number: fnext after finding break " << fnext << endl;

	// Check if the first quark is 1 or 3
	// The factor from the first quark
	fq1=(fnext/( 2factorial(n_g) ))2*factorial(n_g);
	// Special case of 2 quark lines of equal length, then there is no option for q1 (it is 1)
	if( break_after==n_g-break_after ) fq1=0; // 1 is first, no factor from first quark
	//cout << "new_vector_number: fq1 " << fq1 << endl;
	fnext= fnext % (2*factorial(n_g));
	//cout << "new_vector_number: fnext after counting first quark " << fnext << endl;

	// First treat the special case when the ql's switch place
	// This can happen when
	// 1) The 2nd ql gets longer than the first such as
	// emission from 6 in {1,5,2}{3,6,4} -> {3,6,7,4}{1,5,2}
	// 2) The 2nd is one shorter than the first but has the q=1 standing first
	// {3,5,6,2}{1,7,4} ->{1,8,7,4}{3,5,6,2}
	if( ((break_after==n_g-break_after) && place.first==1) or ((break_after==n_g+1-break_after) && fq1>0 && place.first==1)) {
	//cout << "WARNING: this is the hard case" << endl;

	// Locate all partons in old ql
	std::vector< std::pair<int, std::pair<int,int> > > parton_and_place;
	for( int p=1; p<=n_g+2*n_q; p++ ){
	std::pair<int,int> p_place=find_parton(p, oldouble_num, n_q, n_g, n_loop);
	make_pair(p,p_place);
	// note that parton p stands at place p-1
	parton_and_place.push_back( make_pair(p,p_place) ) ;
	}
	//cout << "new_vector_number: reconstructed all parton places" <<endl;

	// To contain new tensor number
	int new_num=0;

	// First calculate the number from the split
	// After the split the number of gluons in the qls are
	// break_after in the 2nd ql
	// This is true for both special cases of form
	// {1,5,2}{3,6,4} -> {3,6,7,4}{1,5,2}
	// {3,5,6,2}{1,7,4} ->{1,8,7,4}{3,5,6,2}
	new_num=(break_after)22*factorial(n_g+1);
	//cout << "new_vector_number: new_num after adding split "<< new_num <<endl;


	// Then calculate the number from the first quark
	// There is a factor 2 from the anti-quark and a factor (n_g+1)! from the gluons
	// if the first factor is 3, which it is in swapings of form
	// {1,5,2}{3,6,4} -> {3,6,7,4}{1,5,2}
	if( (break_after==n_g-break_after) && place.first==1 ) new_num=new_num+2*factorial(n_g+1);
	// In cases of form {3,5,6,2}{1,7,4} ->{1,8,7,4}{3,5,6,2}
	// the factor from the first quark is 0, and
	//cout << "new_vector_number: new_num after first quark "<< new_num <<endl;


	// For the anti-quark qbar=2, check if it stands in first or 2nd ql
	// If it stands in first ql, there is no contribution to the tensor numbers
	// otherwise there is a contribution from all orders of gluons in 2nd ql
	if( parton_and_place.at(2-1).second.first==0) new_num+=factorial(break_after);
	//cout << "new_vector_number: new_num after first anti-quark "<< new_num <<endl;

	// The 2nd anti-quark and 2nd quark gives no contribution

	// Then calculate the contribution to the new tensor number from the gluons
	// Loop over gluons
	for( int p=5; p<=n_g+2*n_q; p++ ){

	// The NEW scalar place of the gluon
	int p_place=parton_and_place.at(p-1).second.second; // contribution from pos in ql
	// If a gluon was standing in the first ql it is afterwards standing in the 2nd
	// All gluons in the new first (old 2nd ql) are therefore standing in front
	// they are n_g-break_after +1 (the +1 is added below)
	if(parton_and_place.at(p-1).second.first==0) p_place+=(n_g-break_after);
	// If the new gluon was inserted before the parton, then the p_place should be increased with one
	// This happens when the parton stands in the new 2nd ql= old first ql (always)
	// or when the gluon is in the new first ql, but after the emitter
	if( parton_and_place.at(p-1).second.first==0 or place.second <=parton_and_place.at(p-1).second.second ) p_place++;// ?? is this right

	//cout << "new_vector_number: looping over gluons, gluon "<< p <<" at new place " << p_place <<endl;


	// See how many smaller stand to the right
	int smaller_right_of_p=0;
	// Check all gluons with smaller number than p to see how many stand to the right
	for( int ps=5; ps<p; ps++ ){
	int ps_place=parton_and_place.at(ps-1).second.second; // Contribution from pos in ql
	if(parton_and_place.at(ps-1).second.first==0) ps_place+=(n_g-break_after); // From first ql
	if( parton_and_place.at(ps-1).second.first==0 or place.second <=parton_and_place.at(ps-1).second.second ) ps_place++; // From new g
	// If the smaller parton ps stand to the right of the parton p
	if( ps_place>p_place ) smaller_right_of_p++;
	}
	//cout << "new_vector_number: Found number of smaller to the right "<< smaller_right_of_p <<endl;

	// Increase the new tensor number for partons to the right
	if( smaller_right_of_p>0 ){
	// If the gluon stands in the first ql (former 2nd), there is also a factor of 2 from the
	// choice of anti-quarks
	// All smaller_right_of_p ways of replacing p with a smaller number contributes a factor
	// factorial(gluons to the right), and there are smaller_right_of_p options
	if(parton_and_place.at(p-1).second.first==1) new_num+=2smaller_right_of_pfactorial(n_g+1-p_place);
	else new_num+=smaller_right_of_p*factorial(n_g+1-p_place);
	}
	//cout << "new_vector_number: new_num after adding info from gluon " << p << " at place " << p_place<< " " << new_num <<endl;

	} // End of loop over gluons

	// Add contribution from the new gluon. The gluon number is larger than every other gluon number,
	// so all gluons standing to the right are smaller
	// In the new first ql the new gluon is at place place.second
	// and there are in total n_g-break_after+1 gluons, so the total number of gluons to the right is
	// break_after+(n_g-break_after+1) -place.second
	int n_right_new=(n_g+1)-place.second;
	// The overall contribution also has a factor 2 from the qbar choice
	// The factor 2 is from the qbar, factorial(n_right_new) from ordering gluons
	new_num+=2n_right_newfactorial(n_right_new);
	//cout << "new_vector_number: new_num after adding info from new gluon " << new_num <<endl;


	return new_num;
	} // End of hard case
	}


	if ( n_q==0 ) place2=find_parton(2, oldouble_num, n_q, n_g, n_loop).second;
	//cout << "The place of 2 was " << place2 << endl;
	//cout.flush();
	// Initially we had q=1, g1,....gn, qbar=2
	// After the insertion at place place we have q=1, g1, ...g(place-1), g_new, g(place), ...gn, qbar=2
	// The number of the new vector can be seen as a sum of 3 factors:
	// f1 = the contribution from q=1, g1, ...g(place-1)
	// fnew= the contribution from the newly inserted gluon
	// f2= the contribution from g(place), ...gn, qbar=2 (this contribution is the same as before the insertion)
	// Of these contributions fnew is easiest to calculate it is simply:
	int fnew=0;
	// If one ql
	if( n_q==1 ) fnew=(n_g - place.second + 1)*factorial(n_g - place.second + 1);
	// Special case of place after last gluon
	else if( n_q==0 && place.second == n_g+1 ) fnew=0;
	// If one closed ql and the parton is to the right of 2
	else if( n_q==0 && place.second > place2 && place.second < n_g+1) fnew=(n_g - place.second )*factorial(n_g - place.second );
	// If the parton is before 2
	else if( n_q==0 && place2 >= place.second ) fnew=(n_g - place.second )*(factorial(n_g - place.second )/2); //Rel order 2 3 fixed


	// as all indices standing to the right are smaller

	//cout << "fnew " << fnew << endl;
	// i.e. the (number of partons with smaller number to the left)(number of partons to the left)!

	// The factor f1 is trickier to obtain. Its initial value (before emission) is
	// sum_{i=1}^(place-1) n_sri (n_p-i)!
	// where n_sri denotes the number of smaller partons standing to the rigt at place i
	// these factors are not known, but can be calculated.
	// Once these factors are known the final value of f1 can also be obtained from
	// sum_{i=1}^(place-1) n_sri (n_p-i+1)!
	// where the extra 1 in the factorial is present as after the emission there is
	// one more parton standing to the right.

	int f1new=0; // To contain the new f1
	int f1old=0; // To contain the old f1

	int scalar_place=0; // The scalar place
	if( place.first!=0 ) scalar_place=break_after;
	scalar_place+=place.second;
	//cout << "Scalar place " << scalar_place << endl;

	if( n_q==2 ) {
	fnew=(n_g + 1- scalar_place )*factorial(n_g + 1- scalar_place );

	// If g inserted in first ql, there is an extra factor of 2 from interchanging quarks
	//if( break_after>=scalar_place-1 ) fnew*=2;
	if( place.first==0 ) fnew*=2;

	//cout << "new_vector_number: scalar_place " << scalar_place << endl;
	//cout << "new_vector_number: fnew " << fnew << endl;
	} // end of if( n_q==2 )



	// The last relevant parton for calculating the contribution to f1
	int last_rel=scalar_place-1;
	// If it's a quark_line and the new gluon is to be inserted after the last gluon
	// then the last place with a q, is not important
	//if(scalar_place==n_g+1) last_rel--;
	// If the g is on the first ql
	//if(place.second==break_after+1) last_rel--;
	// if the g is in the second ql
	//if( place.second+break_after==n_g+1 ) last_rel--;???
	//cout << "last_rel " << last_rel << endl;
	//cout.flush();


	for (int i=1; i<=last_rel; i++){
	//cout << "i " << i << endl;
	// number of partons to the right
	int n_r=n_g-i;
	if( n_q==0 ) n_r--;

	// nsri is the number of smaller partons to the right
	int nsri=0;
	if (n_q==1 ) {
	// Note integer division, how many times do we find the factorial
	nsri=fnext/(factorial(n_r));
	// The contribution to the old f1
	f1old+=nsri*factorial(n_r);
	// The next effective f, to use for finding next nsri
	fnext=fnext%(factorial(n_r));
	// The contribution to the new f1
	f1new+=nsri*(factorial(n_r+1));
	}
	else if (n_q==2) {
	int q_part=1;
	// Keeps track of a possible multiplicative factor from first qbar option
	if ( i <= break_after) q_part=n_q;
	//if (place.first==0) q_part=2;


	// Note integer division, how many times do we find the factorial
	nsri=fnext/(q_part*factorial(n_r));
	// The contribution to the old f1
	f1old+=nsriq_partfactorial(n_r);
	// The next effective f, to use for finding next nsri
	fnext=fnext%(q_part*factorial(n_r));
	// The contribution to the new f1
	f1new+=nsri(q_partfactorial(n_r+1)); //2*2!
	//cout << "i " << i << " nsri " << nsri << endl;

	// Take care of factor from first qbar "when we reach it"
	// if( i == break_after && place.second!=break_after+1 && break_after=n_g-break_after) {
	if( i == break_after ) {
	// Note integer division, how many times do we find the factorial
	nsri=fnext/(factorial(n_r));
	int f1old_before_qbar=f1old;
	// The contribution to the old f1
	f1old+=nsri*factorial(n_r);
	// The next effective f, to use for finding next nsri
	fnext=fnext%(factorial(n_r));
	// if the gluon is inserted just before the qbar the factor from
	// the qbar should be taken care of in the f2 instead
	// (f1old still has to be recalculated as it is used for finding f2 which is unchanged)
	// just added 11 02 03
	//if( !break_after==n_g-break_after ){
	// If the g is NOT inserted just before the qbar
	// i.e. if it is either inserted in 2nd ql or at other place than break_after+1
	if ( !(scalar_place == break_after+1) or place.first ==1){
	// The contribution to the new f1 from the anti-quark
	//f1new+=nsri*(factorial(n_r+1));
	// There are (n_r+1)! options to have 2 before 4
	// and we should only have a contribution if 4 is before 2
	f1new+=nsri*(factorial(n_r+1));
	}
	else { // If the gluon is inserted just before the qbar
	// The factor should be accounted for in fnew instead
	// If the second qbar was 4, this should be compensated for in fnew
	if( f1old_before_qbar!=f1old && place.first ==0){

	// The factor is the number of ways of ordering the gluons in 2nd ql
	// The q and qbar are fixed and do not contribute
	// There is no factor if no gluons to the right ??
	// There are n_r gluons in the second ql (also after emission)
	//if (n_r!=0) Why? even with 0 gluons to the right there is one lower option from qbar
	fnew=fnew+factorial(n_r);
	}
	}
	}

	} // end if (n_q==2)

	// If i is after the place of 2 we have the standard case
	else if(n_q==0 && i >= place2) {
	//cout << "i " << i << " after= 2 " << " n_r " << n_r << endl;

	nsri=fnext/(factorial(n_r));
	f1old+=nsri*factorial(n_r);
	fnext=fnext%(factorial(n_r));
	f1new+=nsri*(factorial(n_r+1));
	}
	// If i is before the place of 2 ordering of 2 and 3 irrelevant
	else if(n_q==0 && i < place2) {
	//cout << "i " << i << " before 2 " << " n_r " << n_r <<endl;

	nsri=fnext/(factorial(n_r)/2);
	f1old+=nsri*(factorial(n_r)/2);
	fnext=fnext%(factorial(n_r)/2);
	f1new+=nsri*(factorial(n_r+1)/2);
	}

	/*
	cout << "f1new " << f1new << endl;
	cout << "f1old " << f1old << endl;
	cout << "fnext " << fnext << endl;
	cout << "f1new " << f1new << endl;
	*/
	}
	/*
	cout << "f1old " << f1old << endl;
	cout << "f1new " << f1new << endl;
	cout << "fnew " << fnew << endl;
	*/

	// Once f1 is known then f2 is simply,
	int f2=0;
	if( n_q==2 ) {
	// For 2 qqbar pairs we have to take care of factor from split
	f2=oldouble_num-f1old-fsplit-fq1;
	}
	else f2=oldouble_num-f1old;
	//cout << "f2 " << f2 << endl;


	// and the contribution from f2 is the same as before (irrespectively of position relative 2)

	// In the case of 2 quarks, if there was a contribution to the tensor number from the split
	// this has to be added back
	// break_after = oldouble_num / (2*factorial(n_g));
	// The contribution was (n_g-break_after)(2factorial(n_g))
	// The special case of ql swapping should already have been taken care of
	if( n_q==2 ) {

	// Calculating new factor from the split

	int fsplitnew=0;
	//fsplitnew=fsplit*(n_g+1);

	// New factor from split after insertion in first ql
	if(place.first==0) {
	// The length of 2nd ql is then still (n_g-break_after)
	// and for each break there are 2*factorial(n_g+1) options
	// The lengths of first and 2nd ql are never equal as 1st was at least as long as 2nd
	// -> factor 2 from option of first q
	fsplitnew=(n_g-break_after)22*factorial(n_g+1);
	}

	// If the insertion was in the second ql (of new length n_g-break_after+1)
	if(place.first==1) {
	// The length of the 2nd ql is now (n_g+1-break_after)
	// all gluon orders matter, and the qbar matters
	// Also the factor fomr the split, coming from all "lower splits"
	// has a factor 2 as the q matters for the lower splits
	fsplitnew=(n_g+1-break_after)22factorial(n_g+1); // 12*2!
	// Unless the lengths are equal there is also a factor from the q
	//if( break_after!=n_g-break_after+1 ) fsplitnew*=2;

	// This is not right one has to loop over splitings some may be eual and some not....
	}


	// New factor from first q in the normal case
	//cout << "fq1 " << fq1 << endl;

	int fq1new=fq1*(n_g+1);
	//cout << "fq1new " << fq1new << endl;

	// If after emission the 2nd ql has the same length as the first
	if(place.first==1 && break_after==n_g+1-break_after){
	fq1new=fq1new/2; // How can this be right?? There should be no factor from first quark then
	}


	f1new=fsplitnew+fq1new+f1new;
	/*
	cout << "fsplitnew " << fsplitnew << endl;
	cout << "fq1new " << fq1new << endl;
	cout << "f1new " << f1new << endl;
	cout << "fnew " << fnew << endl;
	*/
	} // End of n_q==2

	// The total vector number is just the sum of the parts
	//cout << "total " << f1new+fnew+f2 << endl;

	return f1new+fnew+f2;
	}


	// Function for finding the new vector numbers in the basis for n_p+1 partons
	// after radiating a new parton from the parton emitter.
	// The old vector has number oldouble_num, and there were, before emission
	// n_q quarks (+ n_q anti-quarks) and n_g gluons, i.e. n_p=2 n_q+ n_g.
	// For emission from q or qbar there is only one resulting color structure,
	// and -1 is returned in the place of the absent color structure.
	// The second vector, where the new parton is inserted before the emitter
	// comes with a minus sign in the new total amplitude.
	// Function so far only intended for tree-level processes with at most 2 qqbar pairs.
	std::pair<int, int> Col_functions::new_vector_numbers( int oldouble_num, int emitter, int n_q, int n_g, int n_loop ) const{
	//cout << "new_vector_numbers: emission from " << emitter <<" for n_g=" << n_g << " and n_q=" << n_q<< endl;


	if( !(n_q==0 or n_q==1 or n_q==2) or n_loop!=0 ){
	std::cerr << "new_vector_numbers: Function only intended for special case of 0-2 q qbar pair at tree level. For the general case use the general version." << std::endl;
	throw -1;
	}

	// Find the place of the parton
	std::pair<int,int> parton_place=find_parton( emitter, oldouble_num, n_q, n_g, n_loop);
	//cout << "new_vector_numbers: Parton " << emitter << " found at place (" << parton_place.first << " " << parton_place.second << ")" << endl;

	// If the q or qbar is radiating, there is only one term
	if(parton_place.second==0 && n_q>=1)// If the q is radiating
	return std::make_pair(new_vector_number(oldouble_num, std::make_pair(parton_place.first,parton_place.second+1), n_q, n_g, n_loop),-1);
	else if(parton_place.second==n_g+1 && n_q==1)
	return std::make_pair(-1,new_vector_number(oldouble_num, parton_place, n_q, n_g, n_loop));


	// If the emitter is a g there are two different terms
	std::pair<int,int> first_insertion_place, second_insertion_place;
	if( n_q==0 && parton_place.second!=0) {
	first_insertion_place=std::make_pair(parton_place.first, parton_place.second + 1);
	second_insertion_place=std::make_pair(parton_place.first, parton_place.second );
	}
	// Special case that the emitter is standing first
	else if(n_q==0 && parton_place.second==0){
	first_insertion_place=std::make_pair(0, 1);
	second_insertion_place=std::make_pair(0, n_g );
	}
	else if(n_q==1){ // special cases already excluded
	first_insertion_place=std::make_pair(parton_place.first, parton_place.second+1);
	second_insertion_place=std::make_pair(parton_place.first, parton_place.second);
	}
	else if(n_q==2){ // special case of quark already excluded

	// If the emitter is an anti-quark
	if(emitter==2 or emitter==4){
	return std::make_pair(-1,new_vector_number(oldouble_num, parton_place, n_q, n_g, n_loop));
	}

	// Emitter is a gluon, and is thus not standing in the ends
	first_insertion_place=std::make_pair(parton_place.first, parton_place.second+1);

	second_insertion_place=std::make_pair(parton_place.first, parton_place.second);
	}

	int first_tens = new_vector_number(oldouble_num, first_insertion_place, n_q, n_g, n_loop);
	int second_tens = new_vector_number(oldouble_num, second_insertion_place, n_q, n_g, n_loop);

	//cout << "new_vector_numbers: First tensor was " << first_tens << endl;


	return std::make_pair( first_tens, second_tens );
	}

	// Function for finding the new vector numbers in the basis for n_p+1 partons
	// after radiating a new parton from the parton emitter.
	// The old vector is Cs, and after emission it is decomposd in the new basis
	// NewBasis.
	// Returned is the number of the new basis vectors.
	// For emission from q or qbar there is only one resulting color structure,
	// and -1 is returned in the place of the absent color structure.
	// The second vector, where the new parton is inserted before the emitter
	// comes with a minus sign in the new total amplitude.
	std::pair<int, int> Col_functions::new_vector_numbers( const Col_str & Cs, int emitter, const Col_amp & New_basis ) const{

	// Number of new gluon
	int new_g=Cs.n_gluon()+2*Cs.n_quark()+1;

	// To contain the numbers of the new non-zero vectors
	int plus_comp=-1;
	int minus_comp=-1;

	// The vector after emission
	Col_amp Ca=emit_gluon(Cs, emitter, new_g);
	//cout << "new_vector_numbers: Ca " << Ca << endl;
	Poly_vec vec=decompose(Ca, New_basis);

	// Loop over vector entries to identify the non-zero ones
	for (uint veci= 0; veci < New_basis.size(); veci++) {
	if( double_num(vec.at(veci))!=0 ){
	if(double_num(vec.at(veci)) > 0){plus_comp=veci;}
	if(double_num(vec.at(veci)) < 0){minus_comp=veci;}
	}
	}

	std::pair<int, int> res= std::make_pair(plus_comp,minus_comp);

	return res;
	}

	}

File Metadata

Mime Type: text/x-c
Expires: Sat, May 3, 6:45 AM (1 d, 8 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4983128
Default Alt Text: Col_functions.cc (176 KB)

Col_functions.ccNo OneTemporaryActions

Col_functions.ccView Options

File Metadata

Event Timeline

Col_functions.cc
No OneTemporary
Actions

Col_functions.cc
View Options