//-------------------------------------------------------
// example31.C
// X -> p pi- decay
// all values are in GeV
//-------------------------------------------------------
void example21()
{
  double const deg = 180 / 3.1415926;

  // given 4-vectors of decay products
  TLorentzVector proton, pion, X;
  proton.SetXYZM(-0.49, -0.20, 2.11, 0.938);
  pion.SetXYZM(-0.26, -0.05, 0.47, 0.140);

  // setup mother 4-vector
  X = proton + pion;
  cout << "Mass of X = " << X.M() << endl;
  cout << "Pt of X   = " << X.Pt() << endl;
  cout << "Eta of X  = " << X.Eta() << endl;
  cout << "Phi of X  = " << X.Phi() << endl;

  // angle
  double angle = proton.Vect().Angle( pion.Vect() );  // lab frame
  cout << "Angle between proton and pion in lab = " << angle * deg << endl;
  cout << "KE of proton in lab = " << proton.E()-proton.M() << endl;
  cout << "KE of pion   in lab = " << pion.E()-pion.M() << endl;

  // go to CM frame of X
  double betax = X.Px()/X.E();
  double betay = X.Py()/X.E();
  double betaz = X.Pz()/X.E();
  // Boost() function modifies the original 4-vector components
  proton.Boost(-betax, -betay, -betaz);
  pion.Boost(-betax, -betay, -betaz);
  double angle_cm = proton.Vect().Angle( pion.Vect() ); // cm frame

  cout << "Angle between proton and pion in CM  = " << angle_cm * deg << endl;
  cout << "KE of proton in lab = " << proton.E()-proton.M() << endl;
  cout << "KE of pion   in lab = " << pion.E()-pion.M() << endl;
}