MultivariateAnalysis/interface/TMBLorentzVector.h

#ifndef INCLUDE_TMBLORENTZVECTOR
#define INCLUDE_TMBLORENTZVECTOR


// See copyright statements at end of file.

#ifndef INCLUDE_TMBVECTOR3
#include "LJMet/MultivariateAnalysis/interface/TMBVector3.h"
#endif
#ifndef ROOT_TRotation
#include "TRotation.h"
#endif

#ifndef ROOT_TLorentzVector
#include "TLorentzVector.h"
#endif

class TLorentzRotation;

/** 
// TMBLorentzVector is based on ROOT's TLorentzVector, implementing all of its
// interfaces. Changes compared to TLorentzVector:
// * derives from TMBVector3 (instead of having a TVector3 member)
// * mass is stored as member, not energy
// * constructor allows several parameter sets, not just xyzE
// * Rapidity() is rapidity, not pseudo-rapidity
// * Some protections against divisions by zero
// * some virtual methods
//
// NOTE on mass: 
//   by setting an energy < |momentum|, Mag2() is negative.
//   In that case, M() will return -sqrt(-Mag2()).
// \ingroup  reco
*/

class TMBLorentzVector: public TMBVector3 {

public:

    /// Safe indexing of the coordinates when using with matrices, arrays, etc.
   enum ECoordinates { 
       kX = 0, 
       kY = 1, 
       kZ = 2, 
       kE = 3, 
       kNUM_COORDINATES = 4, 
       kSIZE = kNUM_COORDINATES 
   };

    /// Initializer sets for the constructor
   enum EInitializer { 
       kXYZE, ///< the four values are given in Cartesian coords (XYZ) and energy
       kXYZM, ///< the four values are given in Cartesian coords (XYZ) and mass
       kPEtaPhiE, ///< the four values are momentum, eta, phi, and energy
       kPEtaPhiM, ///< the four values are momentum, eta, phi, and mass 
       kPtEtaPhiE, ///< the four values are transverse momentum (see Pt()), eta, phi, and energy
       kPtEtaPhiM ///< the four values are transverse momentum (see Pt()), eta, phi, and mass
   };

   TMBLorentzVector(Double_t a = 0.0, Double_t b = 0.0,
                    Double_t c = 0.0, Double_t d = 0.0, 
                    EInitializer init = kXYZE) 
       : TMBVector3 ()
   {
       /// Initialize a TMBLorentzVector.
       /// The meaning of a,b,c,d depends on the value of init,
       /// see EInitializer
       Set(a,b,c,d,init);
   }

   TMBLorentzVector(const TVector3& vect, Double_t v4, EInitializer init=kXYZE):
      TMBVector3(vect) {
       /// Initialize a TMBLorentzVector with a 3 vector and energy or mass,
       /// depending on the value of init (e.g. it's the mass for init==kPtEtaPhiM)
      switch (init) {
         case kXYZM: case kPEtaPhiM: case kPtEtaPhiM: fM=v4;
         default: SetE(v4);
      }
   }
   TMBLorentzVector(const Double_t * a, EInitializer init=kXYZE) 
    : TMBVector3()
   {
       /// Initialize a TMBLorentzVector with an array (array size not checked!)
       /// The meaning of the array values depends on the value of init,
       /// see EInitializer
       Set(a[0],a[1],a[2],a[3],init);
   }
   TMBLorentzVector(const Float_t * a, EInitializer init=kXYZE) {
       /// Initialize a TMBLorentzVector with an array (array size not checked!)
       /// The meaning of the array values depends on the value of init,
       /// see EInitializer
       Set(a[0],a[1],a[2],a[3],init);
   }

   TMBLorentzVector(const TMBLorentzVector& lv): 
     TMBVector3(lv.Vect()), fM(lv.fM) {
       /// Copy constructor.
     }

   TMBLorentzVector(const TLorentzVector& lv): TMBVector3(lv.Vect()) {
       /// Conversion constructor from a TLorentzVector
      SetE(lv.E()); }

   virtual ~TMBLorentzVector() {}
    /// Destructor


    /// @name Getters
    //@{

   Double_t M() const { 
      /// get mass, see remark on mass in TMBLorentzVector class description
      return fM; }

   virtual Double_t E()  const { /// get energy
      return TMath::Sqrt(fM*fM+Mag32()); }
   Double_t Energy() const { /// get energy
      return E(); }
   Double_t T() const { /// get energy / time component
      return E(); }

   virtual const TMBVector3& Vect() const { /// Get spatial component (const).
      return *this; }
   virtual TMBVector3& Vect() { /// Get spatial component (non-const).
      return *this; }

   void GetXYZT(Double_t *carray) const { 
      /// Getters into an arry (no checking of array bounds!)
      GetXYZ(carray); carray[3]=E(); }
   void GetXYZT(Float_t *carray) const {
      /// Getters into an arry (no checking of array bounds!)
      GetXYZ(carray); carray[3]=E(); }

   Double_t operator () (int i) const;
   Double_t operator [] (int i) const { 
      /// Get components by index (as if it was an array).
      /// You can use ECoordinates to index a dimension.
      return operator ()(i); } 

    //@}

    /// @name Setters
    //@{

   void Set(Double_t a = 0.0, Double_t b = 0.0,
       Double_t c = 0.0, Double_t d = 0.0, 
       EInitializer init = kXYZE) {
       /// Initialize a TMBLorentzVector.
       /// The meaning of a,b,c,d depends on the value of init,
       /// see EInitializer
       switch (init) {
       case kXYZM: SetXYZM(a,b,c,d); break;
       case kPEtaPhiE: SetPEtaPhiE(a,b,c,d); break;
       case kPEtaPhiM: SetPEtaPhiM(a,b,c,d); break;
       case kPtEtaPhiE: SetPtEtaPhiE(a,b,c,d); break;
       case kPtEtaPhiM: SetPtEtaPhiM(a,b,c,d); break;
       default: SetXYZT(a,b,c,d);
       }
   }

   void SetM(Double_t a) { /// set mass
      fM=a; }

   void SetT(Double_t a) { /// set E component. See remark on mass in TMBLorentzVector class description
      Double_t m2=a*a-Mag32(); 
      fM=m2>0?TMath::Sqrt(m2) : -TMath::Sqrt(-m2); 
   }
   void SetE(Double_t a) { /// set E component. See remark on mass in TMBLorentzVector class description
      SetT(a); }

   void SetVect(const TVector3 & vect3) { /// Set spatial component.
      Vect()=vect3; }

   void SetPxPyPzE(Double_t px, Double_t py, Double_t pz, Double_t e) {
      /// set all members using Cartesian coordinates and energy
      SetXYZT(px,py,pz,e); }
   void SetXYZT(Double_t  x, Double_t  y, Double_t  z, Double_t t) {
      /// set all members using Cartesian coordinates and energy
      SetXYZ(x,y,z); SetE(t); }
   void SetXYZM(Double_t  x, Double_t  y, Double_t  z, Double_t m) {
      /// set all members using Cartesian coordinates and mass
      SetXYZ(x,y,z); SetM(m); }
   void SetPtEtaPhiE(Double_t pt, Double_t eta, Double_t phi, Double_t e) {
      /// set all members using transverse momentum, eta, phi and energy
      SetPtEtaPhi(pt,eta,phi); SetE(e); }
   void SetPtEtaPhiM(Double_t pt, Double_t eta, Double_t phi, Double_t m) {
      /// set all members using transverse momentum, eta, phi and mass
      SetPtEtaPhi(pt,eta,phi); fM=m; }
   void SetPEtaPhiE(Double_t p, Double_t eta, Double_t phi, Double_t e) {
      /// set all members using |momentum|, eta, phi and energy
      SetPtEtaPhi(p*TMath::Sin(2.*TMath::ATan(TMath::Exp(-eta))),eta,phi); SetE(e); }
   void SetPEtaPhiM(Double_t p, Double_t eta, Double_t phi, Double_t m) {
      /// set all members using |momentum|, eta, phi and mass
      SetPtEtaPhi(p*TMath::Sin(2.*TMath::ATan(TMath::Exp(-eta))),eta,phi); fM=m; }

    //@}

   TMBLorentzVector & operator = (const TMBLorentzVector & lv) {
      /// assignment
      Vect()=lv.Vect(); fM = lv.fM; return *this; }
   

   TMBLorentzVector   operator +  (const TMBLorentzVector &lv) const {
      /// 4vector addition
      return TMBLorentzVector(Vect()+lv.Vect(), E()+lv.E()); }
   TMBLorentzVector & operator += (const TMBLorentzVector &lv) {
      /// 4vector addition
      Double_t e=E(); Vect()+=lv.Vect(); SetE(e+lv.E()); return *this; }

   TMBLorentzVector   operator -  (const TMBLorentzVector &lv) const {
      /// 4vector substraction
      return TMBLorentzVector(Vect()-lv.Vect(), E()-lv.E()); }
   TMBLorentzVector & operator -= (const TMBLorentzVector &lv) {
      /// 4vector substraction
      Double_t e=E(); Vect()-=lv.Vect(); SetE(e-lv.E()); return *this; }

   TMBLorentzVector operator - () const {
      /// unary minus
      return TMBLorentzVector(-Vect(),-E()); }

   TMBLorentzVector operator * (Double_t a) const {
      /// scaling with real numbers (const)
      return TMBLorentzVector(Vect()*a, E()*a); }
   TMBLorentzVector & operator *= (Double_t a) {
      /// scaling with real numbers (non-const)
      Vect()*=a; fM*=a; return *this; }

   Bool_t operator == (const TMBLorentzVector &lv) const {
      /// Comparison, equality operator. 
      /// Two lorentz vectors are equal if all their componens are equal.
      return (fM==lv.fM && Vect()==lv.Vect()); }
   Bool_t operator != (const TMBLorentzVector &lv) const {
      /// Comparison, in-equality operator. 
      /// Two lorentz vectors are not equal if any of their componens are not equal.
      return !(operator == (lv)); }

   /// Equivalence operator. True if all components are
   /// equivalent within machine precision.
   Bool_t is_equal (const TMBLorentzVector &lv) const;

   Double_t Mag2() const { /// get invariant mass squared. See remark on mass in TMBLorentzVector class description
      return fM*fM; }
   Double_t M2() const { /// get invariant mass squared. See remark on mass in TMBLorentzVector class description
      return Mag2(); }

   Double_t Mag() const { /// get invariant mass. See remark on mass in TMBLorentzVector class description
      return fM;
   }

   Double_t Mt2() const { /// transverse mass squared
      return E()*E() - Z()*Z(); }
   

   Double_t Mt() const { /// Transverse mass. If Mt2() is negative then -sqrt(-Mt2()) is returned.
      Double_t mm=Mt2(); 
      return mm < 0.0 ? -TMath::Sqrt(-mm) : TMath::Sqrt(mm); }

   Double_t Beta() const { 
      /// get beta
      if (E()!=0) return Mag3() / E(); 
      else Warning("Beta()", "E()==0.!"); return 0.; }
   Double_t Gamma() const {
      /// get gamma
      Double_t b = Beta();
      if (b!=1.) return 1.0/TMath::Sqrt(1- b*b);
      else Warning("Gamma()", "Beta()==1.!"); return 0.; }

   Double_t Dot(const TMBLorentzVector &lv) const {
      /// scalar (dot) product
      return E()*lv.E() - Vect().Dot(lv.Vect()); }
   Double_t operator * (const TMBLorentzVector &lv) const {
      /// scalar (dot) product
      return Dot(lv);}

   void SetVectMag(const TVector3 & spatial, Double_t magnitude) {
      /// set vector and mass
      Vect()=spatial; SetE(TMath::Sqrt(magnitude * magnitude + spatial * spatial)); }
   void SetVectM(const TVector3 & spatial, Double_t mass) {
      /// set vector and mass
      SetVectMag(spatial, mass); }
   void SetVectE(const TVector3 & spatial, Double_t e) {
      /// set vector and energy
      Vect()=spatial; SetE(e); }

   Double_t Plus() const { 
      /// Returns t + z.
      /// Related to the positive/negative light-cone component,
      /// which some define this way and others define as (t +/- z)/sqrt(2)
      /// CAVEAT: The values returned are T{+,-}Z. It is known that some authors
      /// find it easier to define these components as (T{+,-}Z)/sqrt(2). Thus
      /// check what definition is used in the physics you're working in and adapt
      /// your code accordingly.
      return E() + Vect().Z(); }
   Double_t Minus() const { 
      /// Returns t - z.
      /// Related to the positive/negative light-cone component,
      /// which some define this way and others define as (t +/- z)/sqrt(2)
      /// CAVEAT: The values returned are T{+,-}Z. It is known that some authors
      /// find it easier to define these components as (T{+,-}Z)/sqrt(2). Thus
      /// check what definition is used in the physics you're working in and adapt
      /// your code accordingly.
      return E() - Vect().Z(); }

   TVector3 BoostVector() const {
      /// Returns the spatial components divided by the time component.
      if (E()!=0.) return 1./E()*Vect();
      else Warning("BoostVector()","E()==0.!"); return TVector3(); }

   void Boost(Double_t x, Double_t y, Double_t z) { 
      /// Lorentz boost by vector xyz
      Boost(TVector3(x,y,z)); }
   void Boost(const TVector3 & b);

   Double_t Rapidity() const { 
      /// get rapidity (i.e. NOT pseudo-rapidity!)
      if (E()!=Pz()) return .5*log( (E()+Pz()) / (E()-Pz()) );
      else Warning("Rapidity", "E()==Pz()!"); return 0.; }


   TMBLorentzVector & operator *= (const TRotation &m) {
      /// rotate using a TRotation
      Vect() *= m; return *this; }
   TMBLorentzVector & Transform(const TRotation &m) {
      /// transform using a TRotation, see TVector3::Transform()
      Vect().Transform(m); return *this; }
   TMBLorentzVector & operator *= (const TLorentzRotation &);
   TMBLorentzVector & Transform(const TLorentzRotation &);
   /// Transformation with HepLorenzRotation.

   operator TLorentzVector() const { 
      /// cast to a TLorentzVector
      return TLorentzVector(X(),Y(),Z(),E()); }

private:
   Double32_t fM; ///< mass
     //ClassDef(TMBLorentzVector,1) // A four vector with (-,-,-,+) metric
};

inline
TMBLorentzVector operator * (Double_t a, const TMBLorentzVector & p) {
    /// global scope scale-by-float operator
   return TMBLorentzVector(a*p.Vect(), a*p.E());
}

/*************************************************************************
* Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
* All rights reserved.                                                  *
*                                                                       *
* For the licensing terms see $ROOTSYS/LICENSE.                         *
* For the list of contributors see $ROOTSYS/README/CREDITS.             *
*************************************************************************/

//------------------------------------------------------------------------------
// Copyright(c) 1995-1997, P.Murat (CDF collaboration, FNAL)
//
// Permission to use, copy, modify and distribute this software and its
// documentation for non-commercial purposes is hereby granted without fee,
// provided that the above copyright notice appears in all copies and
// that both the copyright notice and this permission notice appear in
// the supporting documentation. The authors make no claims about the
// suitability of this software for any purpose.
// It is provided "as is" without express or implied warranty.
// *0001 Mar 29 1999 P.Murat: add forgotten scalar product (dot operator)
//------------------------------------------------------------------------------

#endif // INCLUDE_TMBLORENTZVECTOR
Revision:	1.1
Committed:	Tue Nov 11 23:01:21 2008 UTC (16 years, 5 months ago) by kukartse
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V00-03-01, ZMorph_BASE_20100408, gak040610_morphing, V00-02-02, gak011410, gak010310, ejterm2010_25nov2009, V00-02-01, V00-02-00, gak112409, CMSSW_22X_branch_base, segala101609, V00-01-15, V00-01-14, V00-01-13, V00-01-12, V00-01-11, V00-01-10, gak031009, gak030509, gak022309, gak021209, gak040209, gak012809, V00-01-09, V00-01-08, V00-01-07, V00-01-06, V00-01-05, V00-01-04, V00-00-07, V00-00-06, V00-00-05, V00-00-04, V00-01-03, V00-00-02, V00-00-01, HEAD
Branch point for:	ZMorph-V00-03-01, CMSSW_22X_branch
Log Message:	initial creation of a package for LJMET multivariate analysis
#	User	Rev	Content
1	kukartse	1.1	#ifndef INCLUDE_TMBLORENTZVECTOR
2			#define INCLUDE_TMBLORENTZVECTOR
3
4
5			// See copyright statements at end of file.
6
7			#ifndef INCLUDE_TMBVECTOR3
8			#include "LJMet/MultivariateAnalysis/interface/TMBVector3.h"
9			#endif
10			#ifndef ROOT_TRotation
11			#include "TRotation.h"
12			#endif
13
14			#ifndef ROOT_TLorentzVector
15			#include "TLorentzVector.h"
16			#endif
17
18			class TLorentzRotation;
19
20			/**
21			// TMBLorentzVector is based on ROOT's TLorentzVector, implementing all of its
22			// interfaces. Changes compared to TLorentzVector:
23			// * derives from TMBVector3 (instead of having a TVector3 member)
24			// * mass is stored as member, not energy
25			// * constructor allows several parameter sets, not just xyzE
26			// * Rapidity() is rapidity, not pseudo-rapidity
27			// * Some protections against divisions by zero
28			// * some virtual methods
29			//
30			// NOTE on mass:
31			// by setting an energy < \|momentum\|, Mag2() is negative.
32			// In that case, M() will return -sqrt(-Mag2()).
33			// \ingroup reco
34			*/
35
36			class TMBLorentzVector: public TMBVector3 {
37
38			public:
39
40			/// Safe indexing of the coordinates when using with matrices, arrays, etc.
41			enum ECoordinates {
42			kX = 0,
43			kY = 1,
44			kZ = 2,
45			kE = 3,
46			kNUM_COORDINATES = 4,
47			kSIZE = kNUM_COORDINATES
48			};
49
50			/// Initializer sets for the constructor
51			enum EInitializer {
52			kXYZE, ///< the four values are given in Cartesian coords (XYZ) and energy
53			kXYZM, ///< the four values are given in Cartesian coords (XYZ) and mass
54			kPEtaPhiE, ///< the four values are momentum, eta, phi, and energy
55			kPEtaPhiM, ///< the four values are momentum, eta, phi, and mass
56			kPtEtaPhiE, ///< the four values are transverse momentum (see Pt()), eta, phi, and energy
57			kPtEtaPhiM ///< the four values are transverse momentum (see Pt()), eta, phi, and mass
58			};
59
60			TMBLorentzVector(Double_t a = 0.0, Double_t b = 0.0,
61			Double_t c = 0.0, Double_t d = 0.0,
62			EInitializer init = kXYZE)
63			: TMBVector3 ()
64			{
65			/// Initialize a TMBLorentzVector.
66			/// The meaning of a,b,c,d depends on the value of init,
67			/// see EInitializer
68			Set(a,b,c,d,init);
69			}
70
71			TMBLorentzVector(const TVector3& vect, Double_t v4, EInitializer init=kXYZE):
72			TMBVector3(vect) {
73			/// Initialize a TMBLorentzVector with a 3 vector and energy or mass,
74			/// depending on the value of init (e.g. it's the mass for init==kPtEtaPhiM)
75			switch (init) {
76			case kXYZM: case kPEtaPhiM: case kPtEtaPhiM: fM=v4;
77			default: SetE(v4);
78			}
79			}
80			TMBLorentzVector(const Double_t * a, EInitializer init=kXYZE)
81			: TMBVector3()
82			{
83			/// Initialize a TMBLorentzVector with an array (array size not checked!)
84			/// The meaning of the array values depends on the value of init,
85			/// see EInitializer
86			Set(a[0],a[1],a[2],a[3],init);
87			}
88			TMBLorentzVector(const Float_t * a, EInitializer init=kXYZE) {
89			/// Initialize a TMBLorentzVector with an array (array size not checked!)
90			/// The meaning of the array values depends on the value of init,
91			/// see EInitializer
92			Set(a[0],a[1],a[2],a[3],init);
93			}
94
95			TMBLorentzVector(const TMBLorentzVector& lv):
96			TMBVector3(lv.Vect()), fM(lv.fM) {
97			/// Copy constructor.
98			}
99
100			TMBLorentzVector(const TLorentzVector& lv): TMBVector3(lv.Vect()) {
101			/// Conversion constructor from a TLorentzVector
102			SetE(lv.E()); }
103
104			virtual ~TMBLorentzVector() {}
105			/// Destructor
106
107
108			/// @name Getters
109			//@{
110
111			Double_t M() const {
112			/// get mass, see remark on mass in TMBLorentzVector class description
113			return fM; }
114
115			virtual Double_t E() const { /// get energy
116			return TMath::Sqrt(fM*fM+Mag32()); }
117			Double_t Energy() const { /// get energy
118			return E(); }
119			Double_t T() const { /// get energy / time component
120			return E(); }
121
122			virtual const TMBVector3& Vect() const { /// Get spatial component (const).
123			return *this; }
124			virtual TMBVector3& Vect() { /// Get spatial component (non-const).
125			return *this; }
126
127			void GetXYZT(Double_t *carray) const {
128			/// Getters into an arry (no checking of array bounds!)
129			GetXYZ(carray); carray[3]=E(); }
130			void GetXYZT(Float_t *carray) const {
131			/// Getters into an arry (no checking of array bounds!)
132			GetXYZ(carray); carray[3]=E(); }
133
134			Double_t operator () (int i) const;
135			Double_t operator [] (int i) const {
136			/// Get components by index (as if it was an array).
137			/// You can use ECoordinates to index a dimension.
138			return operator ()(i); }
139
140			//@}
141
142			/// @name Setters
143			//@{
144
145			void Set(Double_t a = 0.0, Double_t b = 0.0,
146			Double_t c = 0.0, Double_t d = 0.0,
147			EInitializer init = kXYZE) {
148			/// Initialize a TMBLorentzVector.
149			/// The meaning of a,b,c,d depends on the value of init,
150			/// see EInitializer
151			switch (init) {
152			case kXYZM: SetXYZM(a,b,c,d); break;
153			case kPEtaPhiE: SetPEtaPhiE(a,b,c,d); break;
154			case kPEtaPhiM: SetPEtaPhiM(a,b,c,d); break;
155			case kPtEtaPhiE: SetPtEtaPhiE(a,b,c,d); break;
156			case kPtEtaPhiM: SetPtEtaPhiM(a,b,c,d); break;
157			default: SetXYZT(a,b,c,d);
158			}
159			}
160
161			void SetM(Double_t a) { /// set mass
162			fM=a; }
163
164			void SetT(Double_t a) { /// set E component. See remark on mass in TMBLorentzVector class description
165			Double_t m2=a*a-Mag32();
166			fM=m2>0?TMath::Sqrt(m2) : -TMath::Sqrt(-m2);
167			}
168			void SetE(Double_t a) { /// set E component. See remark on mass in TMBLorentzVector class description
169			SetT(a); }
170
171			void SetVect(const TVector3 & vect3) { /// Set spatial component.
172			Vect()=vect3; }
173
174			void SetPxPyPzE(Double_t px, Double_t py, Double_t pz, Double_t e) {
175			/// set all members using Cartesian coordinates and energy
176			SetXYZT(px,py,pz,e); }
177			void SetXYZT(Double_t x, Double_t y, Double_t z, Double_t t) {
178			/// set all members using Cartesian coordinates and energy
179			SetXYZ(x,y,z); SetE(t); }
180			void SetXYZM(Double_t x, Double_t y, Double_t z, Double_t m) {
181			/// set all members using Cartesian coordinates and mass
182			SetXYZ(x,y,z); SetM(m); }
183			void SetPtEtaPhiE(Double_t pt, Double_t eta, Double_t phi, Double_t e) {
184			/// set all members using transverse momentum, eta, phi and energy
185			SetPtEtaPhi(pt,eta,phi); SetE(e); }
186			void SetPtEtaPhiM(Double_t pt, Double_t eta, Double_t phi, Double_t m) {
187			/// set all members using transverse momentum, eta, phi and mass
188			SetPtEtaPhi(pt,eta,phi); fM=m; }
189			void SetPEtaPhiE(Double_t p, Double_t eta, Double_t phi, Double_t e) {
190			/// set all members using \|momentum\|, eta, phi and energy
191			SetPtEtaPhi(pTMath::Sin(2.TMath::ATan(TMath::Exp(-eta))),eta,phi); SetE(e); }
192			void SetPEtaPhiM(Double_t p, Double_t eta, Double_t phi, Double_t m) {
193			/// set all members using \|momentum\|, eta, phi and mass
194			SetPtEtaPhi(pTMath::Sin(2.TMath::ATan(TMath::Exp(-eta))),eta,phi); fM=m; }
195
196			//@}
197
198			TMBLorentzVector & operator = (const TMBLorentzVector & lv) {
199			/// assignment
200			Vect()=lv.Vect(); fM = lv.fM; return *this; }
201
202
203			TMBLorentzVector operator + (const TMBLorentzVector &lv) const {
204			/// 4vector addition
205			return TMBLorentzVector(Vect()+lv.Vect(), E()+lv.E()); }
206			TMBLorentzVector & operator += (const TMBLorentzVector &lv) {
207			/// 4vector addition
208			Double_t e=E(); Vect()+=lv.Vect(); SetE(e+lv.E()); return *this; }
209
210			TMBLorentzVector operator - (const TMBLorentzVector &lv) const {
211			/// 4vector substraction
212			return TMBLorentzVector(Vect()-lv.Vect(), E()-lv.E()); }
213			TMBLorentzVector & operator -= (const TMBLorentzVector &lv) {
214			/// 4vector substraction
215			Double_t e=E(); Vect()-=lv.Vect(); SetE(e-lv.E()); return *this; }
216
217			TMBLorentzVector operator - () const {
218			/// unary minus
219			return TMBLorentzVector(-Vect(),-E()); }
220
221			TMBLorentzVector operator * (Double_t a) const {
222			/// scaling with real numbers (const)
223			return TMBLorentzVector(Vect()a, E()a); }
224			TMBLorentzVector & operator *= (Double_t a) {
225			/// scaling with real numbers (non-const)
226			Vect()=a; fM=a; return *this; }
227
228			Bool_t operator == (const TMBLorentzVector &lv) const {
229			/// Comparison, equality operator.
230			/// Two lorentz vectors are equal if all their componens are equal.
231			return (fM==lv.fM && Vect()==lv.Vect()); }
232			Bool_t operator != (const TMBLorentzVector &lv) const {
233			/// Comparison, in-equality operator.
234			/// Two lorentz vectors are not equal if any of their componens are not equal.
235			return !(operator == (lv)); }
236
237			/// Equivalence operator. True if all components are
238			/// equivalent within machine precision.
239			Bool_t is_equal (const TMBLorentzVector &lv) const;
240
241			Double_t Mag2() const { /// get invariant mass squared. See remark on mass in TMBLorentzVector class description
242			return fM*fM; }
243			Double_t M2() const { /// get invariant mass squared. See remark on mass in TMBLorentzVector class description
244			return Mag2(); }
245
246			Double_t Mag() const { /// get invariant mass. See remark on mass in TMBLorentzVector class description
247			return fM;
248			}
249
250			Double_t Mt2() const { /// transverse mass squared
251			return E()E() - Z()Z(); }
252
253
254			Double_t Mt() const { /// Transverse mass. If Mt2() is negative then -sqrt(-Mt2()) is returned.
255			Double_t mm=Mt2();
256			return mm < 0.0 ? -TMath::Sqrt(-mm) : TMath::Sqrt(mm); }
257
258			Double_t Beta() const {
259			/// get beta
260			if (E()!=0) return Mag3() / E();
261			else Warning("Beta()", "E()==0.!"); return 0.; }
262			Double_t Gamma() const {
263			/// get gamma
264			Double_t b = Beta();
265			if (b!=1.) return 1.0/TMath::Sqrt(1- b*b);
266			else Warning("Gamma()", "Beta()==1.!"); return 0.; }
267
268			Double_t Dot(const TMBLorentzVector &lv) const {
269			/// scalar (dot) product
270			return E()*lv.E() - Vect().Dot(lv.Vect()); }
271			Double_t operator * (const TMBLorentzVector &lv) const {
272			/// scalar (dot) product
273			return Dot(lv);}
274
275			void SetVectMag(const TVector3 & spatial, Double_t magnitude) {
276			/// set vector and mass
277			Vect()=spatial; SetE(TMath::Sqrt(magnitude * magnitude + spatial * spatial)); }
278			void SetVectM(const TVector3 & spatial, Double_t mass) {
279			/// set vector and mass
280			SetVectMag(spatial, mass); }
281			void SetVectE(const TVector3 & spatial, Double_t e) {
282			/// set vector and energy
283			Vect()=spatial; SetE(e); }
284
285			Double_t Plus() const {
286			/// Returns t + z.
287			/// Related to the positive/negative light-cone component,
288			/// which some define this way and others define as (t +/- z)/sqrt(2)
289			/// CAVEAT: The values returned are T{+,-}Z. It is known that some authors
290			/// find it easier to define these components as (T{+,-}Z)/sqrt(2). Thus
291			/// check what definition is used in the physics you're working in and adapt
292			/// your code accordingly.
293			return E() + Vect().Z(); }
294			Double_t Minus() const {
295			/// Returns t - z.
296			/// Related to the positive/negative light-cone component,
297			/// which some define this way and others define as (t +/- z)/sqrt(2)
298			/// CAVEAT: The values returned are T{+,-}Z. It is known that some authors
299			/// find it easier to define these components as (T{+,-}Z)/sqrt(2). Thus
300			/// check what definition is used in the physics you're working in and adapt
301			/// your code accordingly.
302			return E() - Vect().Z(); }
303
304			TVector3 BoostVector() const {
305			/// Returns the spatial components divided by the time component.
306			if (E()!=0.) return 1./E()*Vect();
307			else Warning("BoostVector()","E()==0.!"); return TVector3(); }
308
309			void Boost(Double_t x, Double_t y, Double_t z) {
310			/// Lorentz boost by vector xyz
311			Boost(TVector3(x,y,z)); }
312			void Boost(const TVector3 & b);
313
314			Double_t Rapidity() const {
315			/// get rapidity (i.e. NOT pseudo-rapidity!)
316			if (E()!=Pz()) return .5*log( (E()+Pz()) / (E()-Pz()) );
317			else Warning("Rapidity", "E()==Pz()!"); return 0.; }
318
319
320			TMBLorentzVector & operator *= (const TRotation &m) {
321			/// rotate using a TRotation
322			Vect() = m; return this; }
323			TMBLorentzVector & Transform(const TRotation &m) {
324			/// transform using a TRotation, see TVector3::Transform()
325			Vect().Transform(m); return *this; }
326			TMBLorentzVector & operator *= (const TLorentzRotation &);
327			TMBLorentzVector & Transform(const TLorentzRotation &);
328			/// Transformation with HepLorenzRotation.
329
330			operator TLorentzVector() const {
331			/// cast to a TLorentzVector
332			return TLorentzVector(X(),Y(),Z(),E()); }
333
334			private:
335			Double32_t fM; ///< mass
336			//ClassDef(TMBLorentzVector,1) // A four vector with (-,-,-,+) metric
337			};
338
339			inline
340			TMBLorentzVector operator * (Double_t a, const TMBLorentzVector & p) {
341			/// global scope scale-by-float operator
342			return TMBLorentzVector(ap.Vect(), ap.E());
343			}
344
345			/*************************************************************************
346			* Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
347			* All rights reserved. *
348			* *
349			* For the licensing terms see $ROOTSYS/LICENSE. *
350			* For the list of contributors see $ROOTSYS/README/CREDITS. *
351			*************************************************************************/
352
353			//------------------------------------------------------------------------------
354			// Copyright(c) 1995-1997, P.Murat (CDF collaboration, FNAL)
355			//
356			// Permission to use, copy, modify and distribute this software and its
357			// documentation for non-commercial purposes is hereby granted without fee,
358			// provided that the above copyright notice appears in all copies and
359			// that both the copyright notice and this permission notice appear in
360			// the supporting documentation. The authors make no claims about the
361			// suitability of this software for any purpose.
362			// It is provided "as is" without express or implied warranty.
363			// *0001 Mar 29 1999 P.Murat: add forgotten scalar product (dot operator)
364			//------------------------------------------------------------------------------
365
366			#endif // INCLUDE_TMBLORENTZVECTOR