//Example program to use VinciaMHV

#include <iostream>
#include <complex>
#include <cmath>
#include <math.h>
#include "2_Quarks.h"
#include "Pythia8/Basics.h"  //to get 4-vectors from Pythia

using namespace std;
using namespace TwoQuarks;

int main() {

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

////////////////////////create 4 test particles////////////////////////////

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

  //Create vector of input type Vec4 to store 4-momenta, vector of type int for 
  //particle IDs, and vector of type double for helicities
  vector<Vec4> momenta; 
  vector<int> IDs;
  vector<double> helicities;

  Vec4 k0(0., 0., 1., 1.);
  Vec4 k1(0.0, 0.0, -1., 1.);
  Vec4 k2(0.6, 0.8, 0., 1.);
  Vec4 k3(-0.6, -0.8, 0.,  1.);
  
//First particle is quark
  momenta.push_back(k0);                    
  IDs.push_back(1); 
  helicities.push_back(-1.0);
  
//Second particle is gluon
  momenta.push_back(k1);  
  IDs.push_back(21); 
  helicities.push_back(-1.0);
  
//Third particle is quark
  momenta.push_back(k2);  
  IDs.push_back(1); 
  helicities.push_back(-1.0); 

//Fourth particle is gluon
  momenta.push_back(k3);
  IDs.push_back(21);
  helicities.push_back(-1.0);

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

///////////////////////////////Begin output////////////////////////////////

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

//Tell the screen what this program does
  cout << "\nWelcome to the example program.";
  cout << endl;

//display each vector
  cout << "\nDisplay the vectors created for testing purposes. The particles ";
  cout << "have not yet been crossed.\nThe momenta are (px, py, pz, E) (M):\n";
  cout << momenta[0] << momenta[1] << momenta[2] << momenta[3];

  cout << "The IDs are: " << IDs[0] << ", " << IDs[1] << ", ";
  cout << IDs[2] << ", " << IDs[3] << endl;

  cout << "The helicities are: " << helicities[0] << ", ";
  cout << helicities[1] << ", " << helicities[2] << ", ";
  cout << helicities[3] << endl;

//-----------------------------------------------------------------------//
  
 
//Create object of type SPBag(nPartons, momenta, helicities, IDs)
  SPBag mySPBag(4, momenta, helicities, IDs);

//Print the bag of spinor inner products
  cout << "\nPrinting all spinor inner products." << endl;
  mySPBag.printMatrix();
  cout << endl;
  
//Calculate and print full colour amplitude
  cout << "\nThe full colour amplitude, calculated using two steps is: " << endl;
  cout << mySPBag.fullColourAmp(1) << endl;
  
//Create object of type SPBag(nPartons), then initialise
  SPBag mySPBag2(4);
  mySPBag2.init(momenta, helicities, IDs);
  
//Calculate all kinematic colour-ordered amplitudes
  cout << "\nThe kinematic colour-ordered amplitudes are: " << endl;
  mySPBag2.amplitude(1);
  mySPBag2.printAmps();
  
//Calculate and print full colour amplitude
  cout << "\nThe full colour amplitude, calculated using four steps is: " << endl;
  cout << mySPBag2.fullColourAmp() << endl;
  
  return 0;
  
}