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
Error occurred while calculating annotation data.
Log Message:	initial creation of a package for LJMET multivariate analysis
#	Content
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