// lv-3.cc
//              A ZOOM PhysicsVectors tutorial example
//
//              Tutorial example 3 for the LorentzVector class:
//
//              Relativistic Kinematics
// * Determining light-cone nature of a 4-vector
// * Invariant mass calculation
// * Relativistic properties
// * Boosts defined by 4-vectors

// Key lines are marked with                            // <------------------

        // Get the LorentzVector class

#include "ZMutility/ZMenvironment.h"
#include "ZMutility/FixedTypes.h"
#include "PhysicsVectors/LorentzVector.h"
#include "PhysicsVectors/UnitVector.h"
ZM_USING_NAMESPACE( zmpv )
USING( std::cout )
USING( std::endl )

        // This example will use iostream output, so get that --
        // in a portable way.

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

        // This example tests the container methods and as such needs
        // to include the following - I picked vector, any container
        // would be fine.

#include <vector>
using std::vector;

int main() {

// For purposes of this test, it is safe and necessary to ignore ZMexceptions.
// This test exercises some of them but we want to continue beyond each one.
  ZMxPhysicsVectors::setHandler( ZMexIgnoreAlways() );


        // (See tutorial example lv-1 for how to construct and access
        //  LorentzVectors.)

        // generate a few sample vectors

  LorentzVector w1 (  3,  2,  5, Tcomponent (12) );
  LorentzVector w2 ( -1,  3,  4, Tcomponent (5) );
  LorentzVector w3 (  3,  4, -12, Tcomponent (13) );
  LorentzVector w4 ( -3, -1,  4, Tcomponent (5.1) );

// ***** Determining light-cone nature of a 4-vector

        // Here is how to determine if the 4-vector is timelike or spacelike
        //------------------------------------------------------------------

  bool isTL = w1.isTimelike();                              // <-------------

  if (isTL) {
    cout << "Vector w1 = " << w1 << " is Timelike. " << endl;
  }


  bool isSL = w2.isSpacelike();                             // <-------------

  if (isSL) {
    cout << "Whereas vector w2 = " << w2 << " is Spacelike. " << endl;
  }

        // Here is how to determine if the 4-vector is close to lightlike
        //---------------------------------------------------------------

  bool isLL = w3.isLightlike();                             // <-------------
        // returns true if -2* t_**2 *epsilon < restMass2 < 2* t_**2 *epsilon

  if (isLL) {
    cout << "Within epsilon= " << w3.getTolerance() << " w3 = " << w3
         << " is Lightlike." << endl;
  }

        // Also, it's possible to see how far from exact equality restMass2 is:

  cout << w3 << " was " << w3.howLightlike()                // <-------------
       << " away from exactly zero." << endl;

        // Here is how to control tolerance for determining isLightlike
        //-------------------------------------------------------------

  Float8 epsilon = .01;
  isLL = w4.isLightlike( epsilon );                         // <-------------

  if (isLL) {
    cout << "Within epsilon= " << epsilon << " w3 = " << w3
         << " is Lightlike, " << endl << w4.howLightlike()  // <-------------
         << " away from exact equality." << endl;
  }


// ***** Invariant mass calculation

        // Here is how to find the rest mass associated with a LorentzVector
        //------------------------------------------------------------------

  Float8 restmass = w1.restMass();                          // <-------------
  cout << "The rest mass of w1= " << w1 << " is: " << restmass << endl;


        // The rest mass will have the same sign as the t component:

  w1.setT( -w1.t() );
  restmass = w1.restMass();
  cout << "Now, rest mass of w1= " << w1 << " is: " << restmass << endl;

        // If you can use the squared rest mass, that is quicker:

  Float8 restmass2 = w1.restMass2();                        // <-------------
  cout << "Squared rest mass of w1 " << w1 << " is: " << restmass2 << endl;


        // Here is how to find the invariant mass of a pair of 4-vectors
        //--------------------------------------------------------------

        // Note that invariantMass() and invariantMass2() called without
        // any arguments are respectively equivalent to restMass()
        // and restmass2().

  Float8 invariantmass = w2.invariantMass(w3);              // <-------------
  cout << "The invariant mass of " << w2 << " and " << w3
       << " is " << invariantmass << endl;

        // Again, the squared invariant mass is a quicker computation:

  Float8 invariantmass2 = w2.invariantMass2(w3);            // <-------------
  cout << "Squared invariant mass of same vectors is "
       << invariantmass2 << endl;


        // Here is how to find the invariant mass of a collection of 4-vectors
        //--------------------------------------------------------------------

#ifdef DEFECT_STANDARD_CPLUSPLUS        // Prevent any problems with containers
#ifndef NT_MSVCPP

#ifdef NEVER
  vector<LorentzVector> wVec;
     wVec.push_back(w1); wVec.push_back(w2);
     wVec.push_back(w3); wVec.push_back(w4);

  invariantmass  =  invariantMass( wVec );                 // <-------------
  invariantmass2 = invariantMass2( wVec );                 // <-------------

  cout << "The invariant mass of the container of "
       << "4-vectors is " << invariantmass << endl
       << "or, if you prefer, the square of the invariant "
       << "mass is " << invariantmass2<< endl;

#endif
#endif // NT_MSVCPP
#endif // DEFECT_STANDARD_CPLUSPLUS


// ***** Relativistic properties

        // Here is how to find beta and gamma
        //-----------------------------------

  Float8 bet = w1.beta();                                  // <-------------
  cout << "Beta for " << w1 << " is " << bet << endl;

  Float8 gamm = w1.gamma();                                // <-------------
  cout << "Gamma for the same is " << gamm << endl;

        // Here is how to find the pseudorapidity and rapidity
        //----------------------------------------------------

  Float8 pseudoRap = w2.eta();                             // <-------------
  cout << "The pseudo-rapidity eta of " << w2 << " is " << pseudoRap << endl;

  Float8 rapidity = w2.rapidity();                         // <-------------
  cout << "The rapidity of " << w2 << " is " << rapidity << endl;

        // The above two cases both are defined with respect to the
        // z-direction. To find the eta or rapidity with respect to a
        // an arbitrary direction, do the following:

  UnitVector u1( 45, DEGREES, 30, DEGREES);
        // For information on how to construct a UnitVector, see sv-3.cc

  pseudoRap = w2.eta(u1);                                  // <-------------
  cout << "The pseudo-rapidity of " << w2 << " with respect to " << endl
       << u1 << " is " << pseudoRap << endl;

  rapidity = w2.rapidity(u1);                              // <-------------
  cout << "The rapidity of " << w2 << " with respect to " << endl
       << u1 << " is " << rapidity << endl;

        // Finally, just in case you were wondering how to find the rapidity
        // along a vector, we've got that too. To differentiate between
        // that and the above, however, it is known as coLinearRapidity:

  rapidity = w4.coLinearRapidity();                        // <-------------
  cout << "The rapidity along the direction of the vector "
       << w4 << " is " << rapidity << endl;


// ***** Boosts defined by 4-vectors

        // Here is how to find boost to rest frame of a 4-vector
        //------------------------------------------------------

  SpaceVector wBoost = w1.findBoostToCM();                 // <-------------
  cout << "The appropriate boost to take " << w1 << endl
       << " to its rest frame is " << wBoost << endl;

        // w1.findBoostToCM() is equal by definition to w1.findBoostToCM(w1)
        // If w1 were boosted according to wBoost, w1 would then be
        // in its rest frame. There is also a method to give directly this
        // rest frame 4-vector:

  LorentzVector rest4Vec = w1.rest4Vector();               // <-------------
  cout << w1 << "in its rest frame is " << rest4Vec << endl;

        // Along the same lines, there is a method to find the boost
        // vector necessary to boost the rest frame 4-vector
        // back to its original state

  SpaceVector boostV = w1.boostVector();                   // <-------------
  cout << "Boosting the rest frame 4-vector " << rest4Vec
       << " by w1.boostVector() " << endl << boostV << " returns "
       << "the original 4-Vector w1 " << w1 << endl;

 
        // Here is how to find boost to rest frame of two 4-vectors
        //---------------------------------------------------------

  wBoost = w2.findBoostToCM(w3);                           // <-------------
  cout << "The appropriate boost to take " << w2
       << "\nto the rest frame of the center of mass of it and "
       << w3 << "\n is " << wBoost << endl;


        // Here is how to find boost to rest frame of collection of 4-vectors
        //-------------------------------------------------------------------

        // Using the vector<LorentzVectors> declared earlier,

#ifdef DEFECT_STANDARD_CPLUSPLUS        // Prevent any problems with containers
#ifndef NT_MSVCPP

#ifdef NEVER
  boostV = findBoostToCM( wVec );                          // <-------------
  cout << "The appropriate boost for the collection "
       << "of LorentzVectors is " << boostV << endl;

#endif
#endif // NT_MSVCPP
#endif // DEFECT_STANDARD_CPLUSPLUS

  return 0;

} /* end of main */