//This will contain the meat of the functions used to calculate MHV amplitudes

#include "2_Quarks.h"

namespace TwoQuarks {

//------------------------------------------------------------------------//

///////////////////Functions from class SpinorProducts//////////////////////

//------------------------------------------------------------------------//

//Calculate sqrt(k^+) = sqrt(k[0] + k[3])
  complex<double> SpinorProducts::sqrt_k_plus (Vec4 k) {
    kPos = k.pPos();
    if (kPos < 0) {
      return complex<double>(0, sqrt(abs((kPos))));
    }
    else {
      return complex<double>(sqrt(kPos),0);
    }
  }

//Calculate sqrt(k^-) = sqrt(k[0] - k[3])
  complex<double> SpinorProducts::sqrt_k_minus (Vec4 k) {
    kNeg = k.pNeg();
    if (kNeg < 0) {
      return complex<double> (0, sqrt((abs(kNeg))));
    }
    else {
     return complex<double> (sqrt(kNeg), 0);
      
    }
  }
  
//This function calculates exp(i\phi_k) used to calculate the spinor inner product
//exp(i\phi_k) = (k^1 + ik^2)/(\sqrt(k^+))
void SpinorProducts::phasePos (Vec4 k_i, Vec4 k_j) {
  phaseDenom_i = sqrt_k_plus (k_i);
  phaseDenom_j = sqrt_k_plus (k_j);
  if (abs(real(phaseDenom_i)) < TINY && abs(imag(phaseDenom_i)) < TINY) {
    if (abs(real(phaseDenom_j)) < TINY && abs(imag(phaseDenom_j)) < TINY) {
      phase_i = pie;
      phase_j = pie;
    }
    else {
      phase_i = pie;
      real(phase_j) = k_j.px();
      imag(phase_j) = k_j.py();
      phase_j /= phaseDenom_j;
    } 
  }
  else {
    if (abs(real(phaseDenom_j)) < TINY && abs(imag(phaseDenom_j)) < TINY) {
      real(phase_i) = k_i.px();
      imag(phase_i) = k_i.py();
      phase_i /= phaseDenom_i;
      phase_j = pie;
    }
    else {
      real(phase_i) = k_i.px();
      imag(phase_i) = k_i.py();
      phase_i /= phaseDenom_i;
      real(phase_j) = k_j.px();
      imag(phase_j) = k_j.py();
      phase_j /= phaseDenom_j;
    }
  }
}

//This function calculates exp(-i\phi_k) used to calculate the spinor inner product
void SpinorProducts::phaseNeg (Vec4 k_i, Vec4 k_j) {
  phaseDenom_i = sqrt_k_plus (k_i);
  phaseDenom_j = sqrt_k_plus (k_j);
  if (abs(real(phaseDenom_i)) < TINY && abs(imag(phaseDenom_i)) < TINY) {
    if (abs(real(phaseDenom_j)) < TINY && abs(imag(phaseDenom_j)) < TINY) {
      phase_i = pie;
      phase_j = pie;
    }
    else {
      phase_i = pie;
      real(phase_j) = k_j.px();
      imag(phase_j) = -k_j.py();
      phase_j /= phaseDenom_j;
    } 
  }
  else {
    if (abs(real(phaseDenom_j)) < TINY && abs(imag(phaseDenom_j)) < TINY) {
      real(phase_i) = k_i.px();
      imag(phase_i) = -k_i.py();
      phase_i /= phaseDenom_i;
      phase_j = pie;
    }
    else {
      real(phase_i) = k_i.px();
      imag(phase_i) = -k_i.py();
      phase_i /= phaseDenom_i;
      real(phase_j) = k_j.px();
      imag(phase_j) = -k_j.py();
      phase_j /= phaseDenom_j;
    }
  }
}

//This function calculates the spinor inner product <ij> if given 2 Vec4s
//<ij> = sqrt(k_j^+) exp(i\phi_i) - sqrt(k_i^+) exp(i\phi_j)
complex<double> SpinorProducts::spinorProdPos (Vec4 k_i, Vec4 k_j){
  //First calculate phases
  phasePos(k_i, k_j);
  //If k_i^+ == 0
  if (phase_i == pie) { 
    if (phase_j == pie) {
      return 0;
    }
    else {
      return sqrt_k_minus (k_i) * phaseDenom_j;
    }
  }
  else if (phase_j == pie) {
    return - sqrt_k_minus (k_j) * phaseDenom_i;
  }
  else {
    return  phaseDenom_j * phase_i - phaseDenom_i * phase_j;
  }
} 

//This function calculates the spinor inner product <ij> if given 2 Vec4s
//[ij] = sqrt(k_i^+) exp(-i\phi_j) - sqrt(k_j^+) exp(-i\phi_i) 
complex<double> SpinorProducts::spinorProdNeg (Vec4 k_i, Vec4 k_j){
  //First calculate phases
  phaseNeg(k_i, k_j);
  //If k_i^+ == 0
  if (phase_i == pie) { 
    if (phase_j == pie) {
      return 0;
    }
    else {
      return - sqrt_k_minus (k_i) * phaseDenom_j;
    }
  }
  else if (phase_j == pie) {
    return sqrt_k_minus (k_j) * phaseDenom_i;
  }
  else {
    return  phaseDenom_i * phase_j - phaseDenom_j * phase_i;
  }
} 

//------------------------------------------------------------------------//

////////////////////////Functions from class SPBag//////////////////////////

//------------------------------------------------------------------------//

//Combine initialiser with constructor:
  SPBag::SPBag (int nPartonsIn,  vector<Vec4> k_i, vector<double> h_i, vector<int> IDs) {
    Error = false;
    nPartons = nPartonsIn;
    if (nPartons > nPartonsMax) {
      Error = true;
      cout << "Too many particle's. Maximum number of particles is " << nPartonsMax << endl;
    }
    for (i = 0; i < nPartons; i++) {
        spBag[i][i] = 0;
    }
    for (i = 0; i < nAmpsMax; i++) {
      ampBag[i] = 0;
    }
    nAmps = 1;
    nGlue = nPartonsIn - 2;
    
//Create new vectors to store momenta, helicities, IDs,
    vector<Vec4> p_i;
    p_i.reserve(8);
    
//Convert incoming to outgoing
    for (i = 0; i < 2; i++) {
      h_i[i] = -h_i[i];
      if (IDs[i] < 6 && IDs[i] > -6) {
	IDs[i] = -IDs[i];
      }
    }

//Set all outgoing momenta and helicity configurations
    for (i = 0; i < nPartons; i++) {
    //Find quark, put it first
      if (IDs[i] > 0 && IDs[i] < 6) {
	qPos = i;
	ids[0] = IDs[i];
	hel[0] = h_i[i];
	if (i < 2) {
	  p_i.push_back(-k_i[i]);
	}
	else p_i.push_back(k_i[i]);
	break;
      }
    }
    
    
    for (i = 0; i < nPartons; i++) {
    //Find anti-quark, set a variable which knows its place
      if (IDs[i] < 0 && IDs[i] > -6) {
	qBarPos = i;
	break;
      }
    }      
    for (i = 0; i < nPartons; i++) {
    //Find gluons, put them in sequentially
      if (i == qPos || i == qBarPos) {
	continue;
      }
      if (IDs[i] == 21) {
	if (i < qPos && i < qBarPos) {
	  ids[i+1] = IDs[i];
	  hel[i+1] = h_i[i];
	  if (i < 2) {
	  p_i.push_back(-k_i[i]);
	  }
	  else p_i.push_back(k_i[i]);
	}
	else if (i > qPos && i < qBarPos) {
      	  ids[i] = IDs[i];
	  hel[i] = h_i[i];
	  if (i < 2) {
	    p_i.push_back(-k_i[i]);
 	  }
	  else p_i.push_back(k_i[i]);
	}
	else if (i < qPos && i > qBarPos) {
      	  ids[i] = IDs[i];
	  hel[i] = h_i[i];
	  if (i < 2) {
	    p_i.push_back(-k_i[i]);
 	  }
	  else p_i.push_back(k_i[i]);
	}
	//if (i > qBarPos && i > qPos) 
	else { 
	  ids[i-1] = IDs[i];
	  hel[i-1] = h_i[i];
	  p_i.push_back(k_i[i]);
	}
      }
    }
    //Check number of gluons
    if (ids[nGlue] != 21) {
      Error = true;
      cout << "Not enough gluons" << endl;
    }
    
  //Put anti-quark at end of sequence
    ids[nPartons - 1] = IDs[qBarPos];
    hel[nPartons - 1] = h_i[qBarPos];
    if (qBarPos < 2) {
      p_i.push_back(-k_i[qBarPos]);
      
    }
    else p_i.push_back(k_i[qBarPos]);

//If correct, fill ampBag	

//find the number of positive helicity particles
    //else {
      nPos = 0;
      for (n = 0; n < nPartons; n++) { 
        if (hel[n] > 0.0) {
	  nPos += 1;
	}
      }
//If mostly negative helicities, create correct spinor bag
      if (nPos == 2) {
        for (i = 0; i < nPartons; i++) {
	  for (j = i+1; j < nPartons; j++) {
	    spBag[i][j] = spinorProdNeg(p_i[i], p_i[j]);
	    spBag[j][i] = - spBag[i][j];
           }
        }
      }
//If mostly positive helicities, create correct spinor bag 
      else if (nPos == nPartons - 2) {
        for (i = 0; i < nPartons; i++) {
	  for (j = i+1; j < nPartons; j++) {
            spBag[i][j] = spinorProdPos(p_i[i], p_i[j]);
            spBag[j][i] = - spBag[i][j];
	    
	  }
        }
      }        
    

//Else there is an error. Tell user
      else {
	Error = true;
        cout << "\nWrong helicity configuration. Check helicities\n" << endl;
      }
  }


//Initialise the bag
  void SPBag::init (vector<Vec4> k_i, vector<double> h_i, vector<int> IDs) {
    nPos = 0;
    
//Create new vectors to store momenta, helicities, IDs,
    vector<Vec4> p_i;
    p_i.reserve(8);
    
//Convert incoming to outgoing
    for (i = 0; i < 2; i++) {
      h_i[i] = -h_i[i];
      if (IDs[i] < 6 && IDs[i] > -6) {
	IDs[i] = -IDs[i];
      }
    }

//Set all outgoing momenta and helicity configurations
    for (i = 0; i < nPartons; i++) {
    //Find quark, put it first
      if (IDs[i] > 0 && IDs[i] < 6) {
	qPos = i;
	ids[0] = IDs[i];
	hel[0] = h_i[i];
	if (i < 2) {
	  p_i.push_back(-k_i[i]);
	}
	else p_i.push_back(k_i[i]);
	break;
      }
    }
    
    
    for (i = 0; i < nPartons; i++) {
    //Find anti-quark, set a variable which knows its place
      if (IDs[i] < 0 && IDs[i] > -6) {
	qBarPos = i;
	break;
      }
    }      
    for (i = 0; i < nPartons; i++) {
    //Find gluons, put them in sequentially
      if (i == qPos || i == qBarPos) {
	continue;
      }
      if (IDs[i] == 21) {
	if (i < qPos && i < qBarPos) {
	  ids[i+1] = IDs[i];
	  hel[i+1] = h_i[i];
	  if (i < 2) {
	  p_i.push_back(-k_i[i]);
	  }
	  else p_i.push_back(k_i[i]);
	}
	else if (i > qPos && i < qBarPos) {
      	  ids[i] = IDs[i];
	  hel[i] = h_i[i];
	  if (i < 2) {
	    p_i.push_back(-k_i[i]);
 	  }
	  else p_i.push_back(k_i[i]);
	}
	else if (i < qPos && i > qBarPos) {
      	  ids[i] = IDs[i];
	  hel[i] = h_i[i];
	  if (i < 2) {
	    p_i.push_back(-k_i[i]);
 	  }
	  else p_i.push_back(k_i[i]);
	}
	//if (i > qBarPos && i > qPos)
	else {  
	  ids[i-1] = IDs[i];
	  hel[i-1] = h_i[i];
	  p_i.push_back(k_i[i]);
	}
      }
    }
  //Put anti-quark at end of sequence
    ids[nPartons - 1] = IDs[qBarPos];
    hel[nPartons - 1] = h_i[qBarPos];
    if (qBarPos < 2) {
      p_i.push_back(-k_i[qBarPos]); 
    }
    else p_i.push_back(k_i[qBarPos]);

//If correct, fill spBag	

//find the number of positive helicity particles
    //else {
      nPos = 0;
      for (n = 0; n < nPartons; n++) { 
        if (hel[n] > 0.0) {
	  nPos += 1; 
	}
      }
//If mostly negative helicities, create correct spinor bag
      if (nPos == 2) {
        for (i = 0; i < nPartons; i++) {
	  for (j = i+1; j < nPartons; j++) {
	    spBag[i][j] = spinorProdNeg(p_i[i], p_i[j]);
	    spBag[j][i] = - spBag[i][j];
           }
        }
      }
//If mostly positive helicities, create correct spinor bag 
      else if (nPos == nPartons - 2) {
        for (i = 0; i < nPartons; i++) {
	  for (j = i+1; j < nPartons; j++) {
            spBag[i][j] = spinorProdPos(p_i[i], p_i[j]);
            spBag[j][i] = - spBag[i][j];    
	  }
        }
      }        
    

//Else there is an error. Tell user
      else {
	Error = true;
        cout << "\nWrong helicity configuration. Check helicities\n" << endl;
      }
  }
//-----------------------------------------------------------------------//
  void SPBag::amplitude(int type) {
    //If Error then just return error message
    if (Error) {
      cout << "Error" << endl;
      ampBag[0] = pie;
    }
  
  //If good helicity configuration then calculate amplitude
    else if (Error == false) {
    
    //Check type. If type = 0, return single colour-ordered amplitude
      if (type == 0) {
      //Only need 1st entry of ampBag
        nAmps = 1;

      //If mostly positive helicities, return the amplitude
        if (nPos == nPartons - 2) {
        //find the negative gluon
          for (n = 1; n < nPartons - 1; n++) {
            if (hel[n] < 0.0) {
              negGlue = n;
              break;
            }
          }
        //Check whether quark is negative helicity
          if (hel[0] < 0) {
          //Calculate the numerator
            num = pow(spBag[0][negGlue], 3.0) * spBag[nPartons - 1][negGlue];
          //Calculate the denominator.
            //First calculate <qbar q>
            denom = spBag[nPartons - 1][0];
            //Calculate rest of denominator
            for (int k = 0; k < nPartons - 1; k++) {
              denom *= spBag[k][k+1];
            }
          //Calculate amplitude
            if (denom == zero) {
	      ampBag[0] = pie;
	      cout << "\nSingularity" << endl;
	      Error = true;
	    }
	    else {
              amp = num/denom;
              ampBag[0] = amp;
	    }
          }

        //If quark is positive helicity:
          else {
          //Calculate the numerator
            num = pow(spBag[nPartons - 1][negGlue], 3.0) * spBag[0][negGlue];
          //Calculate the denominator.
            //First calculate <qbar q>
            denom = spBag[nPartons - 1][0];
            //Calculate rest of denominator
            for (int k = 0; k < nPartons - 1; k++) {
              denom *= spBag[k][k+1];
            }
          //Calculate amplitude
             if (denom == zero) {
	      ampBag[0] = pie;
	      cout << "\nSingularity" << endl;
	      Error = true;
	    }
	    else {
              amp = num/denom;
              ampBag[0] = amp;
	    }
          }
        }

      //If mostly negative helicities, return the amplitude
        if (nPos == 2) {
        //find the positive gluon
          for (n = 1; n < nPartons - 1; n++) {
            if (hel[n] > 0.0) {
              posGlue = n;
              break;
            }
          }
        //Check whether quark is positive helicity
          if (hel[0] > 0) {
          //Calculate the numerator
            num = pow(spBag[0][posGlue], 3.0) * spBag[nPartons - 1][posGlue];
          //Calculate the denominator.
            //First calculate <qbar q>
            denom = pow(-one, nPartons) * spBag[nPartons - 1][0];
            //Calculate rest of denominator
            for (int k = 0; k < nPartons - 1; k++) {
              denom *= spBag[k][k+1];
            }
          //Calculate amplitude
             if (denom == zero) {
	      ampBag[0] = pie;
	      cout << "\nSingularity" << endl;
	      Error = true;
	    }
	    else{
              amp = num/denom;
              ampBag[0] = amp;
	    }
          }

        //If the quark has negative helicity
          else {
          //Calculate the numerator
            num = pow(spBag[nPartons - 1][posGlue], 3.0) * spBag[0][posGlue];
          //Calculate the denominator.
            //First calculate <qbar q>
            denom = pow(-one, nPartons) * spBag[nPartons - 1][0];
            //Calculate rest of denominator
            for (int k = 0; k < nPartons - 1; k++) {
              denom *= spBag[k][k+1];
            }
          //Calculate amplitude
             if (denom == zero) {
	      ampBag[0] = pie;
	      cout << "\nSingularity" << endl;
	      Error = true;
	    }
	    else {
              amp = num/denom;
              ampBag[0] = amp;
	    }
          }
        }
      }

     //Check type. If type = 1, return array of colour-ordered amplitudes
      else if (type == 1) {
      //We'll use full array ampBag
        nAmps = factorial(nPartons - 2);
        
	//Permute Gluons, i.e. Permute particles 1 to nGlue
	
	//Create list of gluons to permute
	int list[nGlue];
	
	//Put 1, 2, ..., n into list
	for (i = 0; i < nGlue; i++) {
	  list[i] = i + 1;
	}
	
      //Calculate Amplitudes and store in ampBag
	//Create integer for ampBag index
	j = 0;
	
       //If mostly positive helicities, return the amplitude
          if (nPos == nPartons - 2) {
          //find the negative gluon
            for (n = 1; n < nPartons - 1; n++) {
              if (hel[n] < 0.0) {
                negGlue = n;
                break;
              }
            }
          //Check whether quark is negative helicity
            if (hel[0] < 0) {
            //Calculate the numerator
              num = pow(spBag[0][negGlue], 3.0) * spBag[nPartons - 1][negGlue];
	    //Do permutations
              //Calculate the denominator.
              //First calculate <qbar q>
                tempDenom = spBag[nPartons - 1][0];
	      do {
	        //Next calculate <q 1>
	        denom = tempDenom * spBag[0][list[0]];
	        //Next calculate <n qbar>
	        denom *= spBag[list[nGlue-1]][nPartons-1];
                //Calculate rest of denominator
                for (int k = 0; k < nGlue -1 ; k++) {
                  denom *= spBag[list[k]][list[k+1]];
                }
              //Calculate amplitude
                if (denom == zero) {
	          ampBag[j] = pie;
	          cout << "\nSingularity" << endl;
		  Error = true;
	        }
	        else {
                  amp = num/denom;
                  ampBag[j] = amp;
	        }
	        j += 1;
	      }
	       //Continue do loop while permutating numbers in list
	       while (next_permutation(list, list + nGlue) );
	    }

          //If quark is positive helicity:
            else {
            //Calculate the numerator
              num = pow(spBag[nPartons - 1][negGlue], 3.0) * spBag[0][negGlue];
	    //Do permutations
              //Calculate the denominator.
              //First calculate <qbar q>
                tempDenom = spBag[nPartons - 1][0];
	      do {
	        //Next calculate <q 1>
	        denom = tempDenom * spBag[0][list[0]];
	        //Next calculate <n qbar>
	        denom *= spBag[list[nGlue-1]][nPartons-1];
                //Calculate rest of denominator
                for (int k = 0; k < nGlue -1 ; k++) {
                  denom *= spBag[list[k]][list[k+1]];
                }  
	    
              //Calculate amplitude
                if (denom == zero) {
    	        ampBag[j] = pie;
  	          cout << "\nSingularity" << endl;
		  Error = true;
  	        }
	        else {
                  amp = num/denom;
                  ampBag[j] = amp;
	        }
	        j += 1;
	      }
	    //Continue do loop while permutating numbers in list
	      while (next_permutation(list, list + nGlue) );
            }
        //Finish mostly positive helicity section
          }

        //If mostly negative helicities, return the amplitude
          else if (nPos == 2) {
          //find the positive gluon
            for (n = 1; n < nPartons - 1; n++) {
              if (hel[n] > 0.0) {
                posGlue = n;
                break;
              }
            }
          //Check whether quark is positive helicity
            if (hel[0] > 0) {
            //Calculate the numerator
              num = pow(spBag[0][posGlue], 3.0) * spBag[nPartons - 1][posGlue];
	    //Do permutations
              //Calculate the denominator.
                //First calculate <qbar q>
                tempDenom = pow(-one, nPartons) * spBag[nPartons - 1][0];
	      do{
	        //Next calculate <q 1>
	        denom = tempDenom * spBag[0][list[0]];
	        //Next calculate <n qbar>
	        denom *= spBag[list[nGlue-1]][nPartons-1];
                //Calculate rest of denominator
                for (int k = 0; k < nGlue -1 ; k++) {
                  denom *= spBag[list[k]][list[k+1]];
                }  
              
            //Calculate amplitude
               if (denom == zero) {
	        ampBag[j] = pie;
	        cout << "\nSingularity" << endl;
		Error = true;
	        }
	        else{
                  amp = num/denom;
                  ampBag[j] = amp;
	        }
	        j += 1;
	      }
	    //Continue do loop while permutating numbers in list
	      while (next_permutation(list, list + nGlue) );
            }

          //If the quark has negative helicity
            else {
            //Calculate the numerator
              num = pow(spBag[nPartons - 1][posGlue], 3.0) * spBag[0][posGlue];
	    //Do permutations
      	      //Calculate the denominator.
              //First calculate <qbar q>
                tempDenom = pow(-one, nPartons) * spBag[nPartons - 1][0];
	      do{
	        //Next calculate <q 1>
	        denom = tempDenom * spBag[0][list[0]];
	        //Next calculate <n qbar>
	        denom *= spBag[list[nGlue-1]][nPartons-1];
                //Calculate rest of denominator
                for (int k = 0; k < nGlue -1 ; k++) {
                  denom *= spBag[list[k]][list[k+1]];
                }  
              
              //Calculate amplitude
                if (denom == zero) {
	          ampBag[j] = pie;
	          cout << "\nSingularity" << endl;
		  Error = true;
	        }
	        else {
                  amp = num/denom;
                  ampBag[j] = amp;
	        }
	        j += 1;
	      }
	    //Continue do loop while permutating numbers in list
	      while (next_permutation(list, list + nGlue) );
            }
          
      //Finish mostly negative helicity section
        }
	
     //Finish if (type == 1) loop 
      }
      
    //Finish if good helicity structure
    }
    
  //Finish function
  } 
  
//Calculate colour amplitude
  double SPBag::fullColourAmp () {
    //Make colour matrix
    if (Error) {
      cout << "Error made, no amplitude calculated" << endl;
      return -1;
    }
    else if (nPartons == 4) {
 
    //Calculate full colour amplitude
      matrixSquare = 0.;
      for (i = 0; i < 2; i++) {
	matrixSquare += norm(ampBag[i]) * QuarkTwoGlueColour[i][i];
        for (n = 0; n < i; n++) {
  	  matrixSquare += 2 * (real(ampBag[i]) * real(ampBag[n]) + 
	      imag(ampBag[i])*imag(ampBag[n])) * QuarkTwoGlueColour[i][n]; 
        }
      }
      return matrixSquare;
    }
   
    else if (nPartons == 5) {
      matrixSquare = 0.;
      for (i = 0; i < 6; i++) {
	matrixSquare += norm(ampBag[i]) * QuarkThreeGlueColour[i][i];
        for (n = 0; n < i; n++) {
  	  matrixSquare += 2 * (real(ampBag[i]) * real(ampBag[n]) + 
	      imag(ampBag[i])*imag(ampBag[n])) * QuarkThreeGlueColour[i][n]; 
        }
      }
      return matrixSquare;
    }
    
    else if (nPartons == 6) {
      matrixSquare = 0.;
      for (i = 0; i < 24; i++) {
	matrixSquare += norm(ampBag[i]) * QuarkFourGlueColour[i][i];
        for (n = 0; n < i; n++) {
  	  matrixSquare += 2 * (real(ampBag[i]) * real(ampBag[n]) + 
	      imag(ampBag[i])*imag(ampBag[n])) * QuarkFourGlueColour[i][n]; 
        }
      }
      return matrixSquare;
    }
    
    
    else if (nPartons == 7) {
      matrixSquare = 0.;
      for (i = 0; i < 120; i++) {
	matrixSquare += norm(ampBag[i]) * QuarkFiveGlueColour[i][i];
        for (n = 0; n < i; n++) {
	  matrixSquare += 2 *(real(ampBag[i]) * real(ampBag[n]) + 
	      imag(ampBag[i])*imag(ampBag[n])) * QuarkFiveGlueColour[i][n]; 
        }
      }
      return matrixSquare; 
    }
    
    
    else return -1.;
    
      
  }
  
//Overlead colour amp squared
  double SPBag::fullColourAmp (int type) {
     //If Error then just return error message
    if (Error) {
      cout << "Error" << endl;
      ampBag[0] = pie;
    }
  
  //If good helicity configuration then calculate amplitude
    else if (Error == false) {
    
    //Specify number of gluons
      nGlue = nPartons - 2;
      
    //Check type. If type = 0, return single colour-ordered amplitude
      if (type == 0) {
      //Only need 1st entry of ampBag
        nAmps = 1;

      //If mostly positive helicities, return the amplitude
        if (nPos == nPartons - 2) {
        //find the negative gluon
          for (n = 1; n < nPartons - 1; n++) {
            if (hel[n] < 0.0) {
              negGlue = n;
              break;
            }
          }
        //Check whether quark is negative helicity
          if (hel[0] < 0) {
          //Calculate the numerator
            num = pow(spBag[0][negGlue], 3.0) * spBag[nPartons - 1][negGlue];
          //Calculate the denominator.
            //First calculate <qbar q>
            denom = spBag[nPartons - 1][0];
            //Calculate rest of denominator
            for (int k = 0; k < nPartons - 1; k++) {
              denom *= spBag[k][k+1];
            }
          //Calculate amplitude
            if (denom == zero) {
	      ampBag[0] = pie;
	      cout << "\nSingularity" << endl;
	      Error = true;
	    }
	    else {
              amp = num/denom;
              ampBag[0] = amp;
	    }
          }

        //If quark is positive helicity:
          else {
          //Calculate the numerator
            num = pow(spBag[nPartons - 1][negGlue], 3.0) * spBag[0][negGlue];
          //Calculate the denominator.
            //First calculate <qbar q>
            denom = spBag[nPartons - 1][0];
            //Calculate rest of denominator
            for (int k = 0; k < nPartons - 1; k++) {
              denom *= spBag[k][k+1];
            }
          //Calculate amplitude
             if (denom == zero) {
	      ampBag[0] = pie;
	      cout << "\nSingularity" << endl;
	      Error = true;
	    }
	    else {
              amp = num/denom;
              ampBag[0] = amp;
	    }
          }
        }

      //If mostly negative helicities, return the amplitude
        if (nPos == 2) {
        //find the positive gluon
          for (n = 1; n < nPartons - 1; n++) {
            if (hel[n] > 0.0) {
              posGlue = n;
              break;
            }
          }
        //Check whether quark is positive helicity
          if (hel[0] > 0) {
          //Calculate the numerator
            num = pow(spBag[0][posGlue], 3.0) * spBag[nPartons - 1][posGlue];
          //Calculate the denominator.
            //First calculate <qbar q>
            denom = pow(-one, nPartons) * spBag[nPartons - 1][0];
            //Calculate rest of denominator
            for (int k = 0; k < nPartons - 1; k++) {
              denom *= spBag[k][k+1];
            }
          //Calculate amplitude
             if (denom == zero) {
	      ampBag[0] = pie;
	      cout << "\nSingularity" << endl;
	      Error = true;
	    }
	    else{
              amp = num/denom;
              ampBag[0] = amp;
	    }
          }

        //If the quark has negative helicity
          else {
          //Calculate the numerator
            num = pow(spBag[nPartons - 1][posGlue], 3.0) * spBag[0][posGlue];
          //Calculate the denominator.
            //First calculate <qbar q>
            denom = pow(-one, nPartons) * spBag[nPartons - 1][0];
            //Calculate rest of denominator
            for (int k = 0; k < nPartons - 1; k++) {
              denom *= spBag[k][k+1];
            }
          //Calculate amplitude
             if (denom == zero) {
	      ampBag[0] = pie;
	      cout << "\nSingularity" << endl;
	      Error = true;
	    }
	    else {
              amp = num/denom;
              ampBag[0] = amp;
	    }
          }
        }
      //Return the amplitude
        if (Error) {
	  cout << "Error, no amplitude calculated" << endl;
	  return -1;
	}
        else if (nPartons == 4) {
	  return real(conj(ampBag[0]) * ampBag[0]) * LC4; 
	}
	else if (nPartons == 5) {
	  return real(conj(ampBag[0]) * ampBag[0]) * LC5; 
	}
	else if (nPartons == 6) {
	  return real(conj(ampBag[0]) * ampBag[0]) * LC6; 
	}
	else if (nPartons == 7) {
	  return real(conj(ampBag[0]) * ampBag[0]) * LC7; 
	}
      }

     //Check type. If type = 1, return array of colour-ordered amplitudes
      else if (type == 1) {
      //We'll use full array ampBag
        nAmps = factorial(nPartons - 2);
        
	//Permute Gluons, i.e. Permute particles 1 to (nPartons - 1)
	
	//Create list of gluons to permute
	int list[nGlue];
	
	//Put 1, 2, ..., n into list
	for (i = 0; i < nGlue; i++) {
	  list[i] = i + 1;
	}
	
      //Calculate Amplitudes and store in ampBag
	//Create integer for ampBag index
	j = 0;
	
       //If mostly positive helicities, return the amplitude
          if (nPos == nPartons - 2) {
          //find the negative gluon
            for (n = 1; n < nPartons - 1; n++) {
              if (hel[n] < 0.0) {
                negGlue = n;
                break;
              }
            }
          //Check whether quark is negative helicity
            if (hel[0] < 0) {
            //Calculate the numerator
              num = pow(spBag[0][negGlue], 3.0) * spBag[nPartons - 1][negGlue];
	    //Do permutations
              //Calculate the denominator.
              //First calculate <qbar q>
                tempDenom = spBag[nPartons - 1][0];
	      do {
	        //Next calculate <q 1>
	        denom = tempDenom * spBag[0][list[0]];
	        //Next calculate <n qbar>
	        denom *= spBag[list[nGlue-1]][nPartons-1];
                //Calculate rest of denominator
                for (int k = 0; k < nGlue -1 ; k++) {
                  denom *= spBag[list[k]][list[k+1]];
                }
              //Calculate amplitude
                if (denom == zero) {
	          ampBag[j] = pie;
	          cout << "\nSingularity" << endl;
		  Error = true;
	        }
	        else {
                  amp = num/denom;
                  ampBag[j] = amp;
	        }
	        j += 1;
	      }
	       //Continue do loop while permutating numbers in list
	       while (next_permutation(list, list + nGlue) );
	    }

          //If quark is positive helicity:
            else {
            //Calculate the numerator
              num = pow(spBag[nPartons - 1][negGlue], 3.0) * spBag[0][negGlue];
	    //Do permutations
              //Calculate the denominator.
              //First calculate <qbar q>
                tempDenom = spBag[nPartons - 1][0];
	      do {
	        //Next calculate <q 1>
	        denom = tempDenom * spBag[0][list[0]];
	        //Next calculate <n qbar>
	        denom *= spBag[list[nGlue-1]][nPartons-1];
                //Calculate rest of denominator
                for (int k = 0; k < nGlue -1 ; k++) {
                  denom *= spBag[list[k]][list[k+1]];
                }  
	    
              //Calculate amplitude
                if (denom == zero) {
    	        ampBag[j] = pie;
  	          cout << "\nSingularity" << endl;
		  Error = true;
  	        }
	        else {
                  amp = num/denom;
                  ampBag[j] = amp;
	        }
	        j += 1;
	      }
	    //Continue do loop while permutating numbers in list
	      while (next_permutation(list, list + nGlue) );
            }
        //Finish mostly positive helicity section
          }

        //If mostly negative helicities, return the amplitude
          else if (nPos == 2) {
          //find the positive gluon
            for (n = 1; n < nPartons - 1; n++) {
              if (hel[n] > 0.0) {
                posGlue = n;
                break;
              }
            }
          //Check whether quark is positive helicity
            if (hel[0] > 0) {
            //Calculate the numerator
              num = pow(spBag[0][posGlue], 3.0) * spBag[nPartons - 1][posGlue];
	    //Do permutations
              //Calculate the denominator.
                //First calculate <qbar q>
                tempDenom = pow(-one, nPartons) * spBag[nPartons - 1][0];
	      do{
	        //Next calculate <q 1>
	        denom = tempDenom * spBag[0][list[0]];
	        //Next calculate <n qbar>
	        denom *= spBag[list[nGlue-1]][nPartons-1];
                //Calculate rest of denominator
                for (int k = 0; k < nGlue -1 ; k++) {
                  denom *= spBag[list[k]][list[k+1]];
                }  
              
            //Calculate amplitude
               if (denom == zero) {
	        ampBag[j] = pie;
	        cout << "\nSingularity" << endl;
		Error = true;
	        }
	        else{
                  amp = num/denom;
                  ampBag[j] = amp;
	        }
	        j += 1;
	      }
	    //Continue do loop while permutating numbers in list
	      while (next_permutation(list, list + nGlue) );
            }

          //If the quark has negative helicity
            else {
            //Calculate the numerator
              num = pow(spBag[nPartons - 1][posGlue], 3.0) * spBag[0][posGlue];
	    //Do permutations
      	      //Calculate the denominator.
              //First calculate <qbar q>
                tempDenom = pow(-one, nPartons) * spBag[nPartons - 1][0];
	      do{
	        //Next calculate <q 1>
	        denom = tempDenom * spBag[0][list[0]];
	        //Next calculate <n qbar>
	        denom *= spBag[list[nGlue-1]][nPartons-1];
                //Calculate rest of denominator
                for (int k = 0; k < nGlue -1 ; k++) {
                  denom *= spBag[list[k]][list[k+1]];
                }  
              
              //Calculate amplitude
                if (denom == zero) {
	          ampBag[j] = pie;
	          cout << "\nSingularity" << endl;
		  Error = true;
	        }
	        else {
                  amp = num/denom;
                  ampBag[j] = amp;
	        }
	        j += 1;
	      }
	    //Continue do loop while permutating numbers in list
	      while (next_permutation(list, list + nGlue) );
            }
          
      //Finish mostly negative helicity section
        }
        
       //Make colour matrix
          if (Error) {
            cout << "Error made, no amplitude calculated" << endl;
            return -1;
          }
        else if (nPartons == 4) {
 
        //Calculate full colour amplitude
          matrixSquare = 0.;
          for (i = 0; i < 2; i++) {
	    matrixSquare += norm(ampBag[i]) * QuarkTwoGlueColour[i][i]; 
            for (n = 0; n < i; n++) {
  	       matrixSquare += 2 * (real(ampBag[i]) * real(ampBag[n]) + 
	      imag(ampBag[i])*imag(ampBag[n])) * QuarkTwoGlueColour[i][n]; 
            }
          }
          return matrixSquare;
        }
   
        else if (nPartons == 5) {
          matrixSquare = 0.;
          for (i = 0; i < 6; i++) {
	    matrixSquare += norm(ampBag[i]) * QuarkThreeGlueColour[i][i];
            for (n = 0; n < i; n++) {
  	      matrixSquare += 2 * (real(ampBag[i]) * real(ampBag[n]) + 
	      imag(ampBag[i])*imag(ampBag[n])) * QuarkThreeGlueColour[i][n]; 
            }
          }
          return matrixSquare;
        }
    
        else if (nPartons == 6) {
          matrixSquare = 0.;
          for (i = 0; i < 24; i++) {
	    matrixSquare += norm(ampBag[i]) * QuarkFourGlueColour[i][i]; 
            for (n = 0; n < i; n++) {
              matrixSquare += 2 * (real(ampBag[i]) * real(ampBag[n]) + 
	      imag(ampBag[i])*imag(ampBag[n])) * QuarkFourGlueColour[i][n]; 
            }
          }
          return matrixSquare;
        }
        
        else if (nPartons == 7) {
	  matrixSquare = 0.;
	  for (i = 0; i < 120; i++) {
	    matrixSquare += norm(ampBag[i]) * QuarkFiveGlueColour[i][i];
            for (n = 0; n < i; n++) {
              matrixSquare += 2 * (real(ampBag[i]) * real(ampBag[n]) + 
	      imag(ampBag[i])*imag(ampBag[n])) * QuarkFiveGlueColour[i][n]; 
            }
          }
          return matrixSquare;
        }
    
        else return -1.;
   
	
     //Finish if (type == 1) loop 
      }
      
    //Finish if good helicity structure
    }
    
  //Finish function
  }

//Finish namespace
}