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

#include "ZMutility/ZMenvironment.h"

#include "PhysicsVectors/SpaceVector.h"
#include "PhysicsVectors/LorentzVector.h"
ZM_USING_NAMESPACE( zmpv )

#include "ZMutility/cmath"
using CMATH_NAMESPACE::fabs;

#include "ZMutility/iostream"
using std::cout;
using std::cin;
using std::endl;

#include "ZMutility/iomanip"

#include <vector>
using std::vector;

//bool problems (false);
//LorentzVector vans;
//LorentzVector v1 ( 1., 2., 3., Tcomponent(4.) );
//LorentzVector v2 ( Tcomponent(5), 1., 2, 3. );
//const LorentzVector v1c ( Tcomponent(6), 2, 3., 4. );
//SpaceVector vspace (5,7,2);
//LorentzVector v3( vspace, 30);
//LorentzVector v4( 20, vspace);
//LorentzVector v5( vspace, Tcomponent(40) );
//LorentzVector v6( Tcomponent(50), vspace );
//LorentzVector v7( Tcomponent(25) );
//LorentzVector v8( -10 );
//SpaceVector vg = v1.v();
//SpaceVector svans(1,2,3);
//SpaceVector vx ( 3, -4, 12);
//const LorentzVector v3c( 3, 5, 7, Tcomponent(9) );
//const LorentzVector v4c( 6, 2, 8, Tcomponent(7) );
//double toler;

bool close ( double a, double b ) {
  return ( fabs(a-b) < .0000000001 );
}


void boosts();


int main() {

  boosts();

  return 0;

} /* End of lvtest main() */

void boosts() {

  int i;

  cout << "\n --- Boosts\n\n";

  cout << "Single particle boosts" << endl;

  LorentzVector v0(0.0,0.0,0.0,Tcomponent(120.0));
  LorentzVector vpx,vpy,vmx,vmy;
  SpaceVector beta;
  double magbeta;

  for( i=0; i<10; ++i ) {
    magbeta = (double)i/10.0 + 0.05;
    cout << "\nStarting with v = " << v0 << " and magbeta = " << magbeta  << endl;
    vpx  = boostXOf(v0,magbeta);
    cout << " boost along  x to get " << vpx << endl;
    vmx = boostXOf(vpx,-magbeta);
    cout << " boost along -x to get " << vmx << endl;

    cout << "\nStarting with v = " << v0 << " and magbeta = " << magbeta << endl;
    vpx  = boostXOf(v0,magbeta);
    cout << " boost along  x to get " << vpx << endl;
    vpy  = boostYOf(vpx,magbeta);
    cout << " boost along  y to get " << vpy << endl;
    beta = vpy.findBoostToCM();
    cout << " beta for boost back to rest frame is: " << beta << endl;
    vpx = vpy.boost(beta);
    cout << " boost back to rest frame to get " << vpx << endl;

    cout << "\nStarting with v = " << v0 << " and magbeta = " << magbeta << endl;
    vpx  = boostXOf(v0,magbeta);
    cout << " boost along  x to get " << vpx << endl;
    vpy  = boostYOf(vpx,magbeta);
    cout << " boost along  y to get " << vpy << endl;
    vmx = boostXOf(vpy,-magbeta);
    cout << " boost along -x to get " << vmx << endl;
    beta = vmx.findBoostToCM();
    cout << " beta for boost back to rest frame is: " << beta << endl;
    vmy = vmx.boost(beta);
    cout << " boost back to rest frame to get " << vmy << endl;
  }

  cout << "\nTwo particle boosts" << endl;

  LorentzVector va( 1.0, 0.0, 0.0,Tcomponent(120.0));
  LorentzVector vb(-1.0, 0.0, 0.0,Tcomponent(120.0));
  LorentzVector vap,vapp,vappp,vbp,vbpp,vbppp;
  SpaceVector sva,svappp;

  for( i=0; i<10; ++i ) {
    magbeta = (double)i/10.0 + 0.05;
    cout << "\nStarting with va = " << va << ", vb = " << vb << endl;
    cout << " and magbeta = " << magbeta << endl;
    vap = boostXOf(va,magbeta);
    vbp = boostXOf(vb,magbeta);
    cout << " boost along x to get va = " << vap << endl;
    cout << "                  and vb = " << vbp << endl;
    vapp = boostYOf(vap,magbeta);
    vbpp = boostYOf(vbp,magbeta);
    cout << " boost along y to get va = " << vapp << endl;
    cout << "                  and vb = " << vbpp << endl;
    beta = vapp.findBoostToCM(vbpp);
    cout << " beta for boost back to rest frame is: " << beta << endl;
    vappp = vapp.boost(beta);
    vbppp = vbpp.boost(beta);
    cout << " boost back to rest frame to get va = " << va << endl;
    cout << " boost back to rest frame to get vb = " << vb << endl;

    sva = va.v();
    svappp = vappp.v();
    double costheta = sva.dot(svappp);
    cout << " for magbeta = " << magbeta << ", costheta = " << costheta << endl;
  }

   cout << "\nThree particle boosts" << endl;

  LorentzVector vc( 10.0, 10.0,-15.0,Tcomponent(120.0));
  LorentzVector vd(  5.0, 10.0, 25.0,Tcomponent(120.0));
  LorentzVector ve(-15.0,-20.0,-10.0,Tcomponent(120.0));
  LorentzVector vcp,vcpp,vcppp;
  LorentzVector vdp,vdpp,vdppp;
  LorentzVector vep,vepp,veppp;
  SpaceVector cmsum;

  for( i=0; i<10; ++i ) {
    magbeta = (double)i/10.0 + 0.05;
    cout << "\nStarting with vc = " << vc << ", vd = " << vd << endl;
    cout << "                ve = " << ve << ", and magbeta = " << magbeta << endl;
    vcp = boostXOf(vc,magbeta);
    vdp = boostXOf(vd,magbeta);
    vep = boostXOf(ve,magbeta);
    cout << " boost along x to get vc = " << vcp << endl;
    cout << "                  and vd = " << vdp << endl;
    cout << "                  and ve = " << vep << endl;
    vcpp = boostYOf(vcp,magbeta);
    vdpp = boostYOf(vdp,magbeta);
    vepp = boostYOf(vep,magbeta);
    cout << " boost along y to get vc = " << vcpp << endl;
    cout << "                  and vd = " << vdpp << endl;
    cout << "                  and ve = " << vepp << endl;
    beta = -(vcpp.getV()+vdpp.getV()+vepp.getV())/(vcpp.t()+vdpp.t()+vepp.t());
    cout << " beta for boost back to rest frame is: " << beta << endl;
    vcppp = vcpp.boost(beta);
    vdppp = vdpp.boost(beta);
    veppp = vepp.boost(beta);
    cout << " boost back to rest frame to get vc = " << vcppp << endl;
    cout << " boost back to rest frame to get vd = " << vdppp << endl;
    cout << " boost back to rest frame to get ve = " << veppp << endl;
    cmsum = vcppp.getV()+vdppp.getV()+veppp.getV();
    cout << " momentum sum in putative three-body CM is " << cmsum << endl;
  }
  cout << endl;

  cout << "\nMany particle boosts" << endl;

  LorentzVector vm[11];
  double MHiggs  = 120.0;
  double Mproton =   0.93827231;
  double Mpion   =   0.13956995;
  double Mmuon   =   0.105658389;
  double Eproton, Epion, EHiggs, Emuon, pmuon;
  double pmux, pmuz, cosTheta;
  int many;

  cout << "\nGenerate a Tevatron-style pbar-p collision that produces one " << endl;
  cout << " 120 GeV Higgs and a gang of pions." << endl;

  for( int j=0; j<20; ++j ) {
    Eproton = sqrt(Mproton*Mproton+900.0*900.0);
    vm[0].set(0.0, 0.0, 900.0, Tcomponent(Eproton));
    vm[1].set(0.0, 0.0,-900.0, Tcomponent(Eproton));
    Epion = sqrt(Mpion*Mpion+100.0*100.0);
    vm[2].set(100.0,   0.0,   0.0, Tcomponent(Epion));
    vm[3].set(  0.0, 100.0,   0.0, Tcomponent(Epion));
    vm[4].set(  0.0,   0.0, 100.0, Tcomponent(Epion));
    Epion = sqrt(Mpion*Mpion+200.0*200.0);
    vm[5].set(200.0,   0.0,   0.0, Tcomponent(Epion));
    vm[6].set(  0.0, 200.0,   0.0, Tcomponent(Epion));
    vm[7].set(  0.0,   0.0, 200.0, Tcomponent(Epion));
    EHiggs = sqrt(MHiggs*MHiggs+13.0*13.0);
    vm[8].set(  5.0,  12.0,   0.0, Tcomponent(EHiggs));
    if( j == 0 ) {
      cout << " momenta of the colliding beams: " << vm[0] << " and " << vm[1] << endl;
      cout << " momentum of the Higgs:          " << vm[8] << endl;
  }
    many = 9;

// Go to the Higgs rest frame

    if( j == 0 ) {
      cout << "\n Go to the Higgs rest frame and let the Higgs decay into two muons" << endl;
    }
    beta = vm[8].getV()/vm[8].t();
    for( i=0; i<many; ++i ) {
      vm[i].boost(beta);
    }

// Let the Higgs decay into two muons

    Emuon = 0.5*MHiggs;
    pmuon = sqrt(Emuon*Emuon-Mmuon*Mmuon);
    cosTheta = 0.95 - 0.1*(double)j;
    pmuz = pmuon*cosTheta;
    pmux = pmuon*sqrt(1.0-cosTheta*cosTheta);
    vm[9].set( pmux, 0.0, pmuz, Tcomponent(Emuon));
    vm[10].set(-pmux, 0.0,-pmuz, Tcomponent(Emuon));
    many = many + 2;

// and go back to the lab frame

    if( j == 0 ) {
      cout << " then take the whole mess (including the muons) back to the lab frame" << endl;
    }
    for( i=0; i<many; ++i ) {
      vm[i].boost(-beta);
    }

// Check that we really are back in the lab

    if( j == 0 ) {
      cout << "\n momenta of the colliding beams: " << vm[0] << " and " << vm[1] << " GeV" << endl;
    }
    cout << " for cosTheta = " << cosTheta << endl;
    cout << " momenta of the two muons: " << vm[9] << " GeV" << endl;
    cout << "                       and " << vm[10] << " GeV" << endl;
    cout << " effective mass of the two muons is " << vm[9].invariantMass(vm[10]) << " GeV" << endl;
    cout << endl;
  }

} // end of boosts()