BaseParticlePropagator/interface/RawParticle.h

0001 #ifndef CommonTools_BaseParticlePropagator_RawParticle_h
0002 #define CommonTools_BaseParticlePropagator_RawParticle_h
0003
0004 #include "DataFormats/Math/interface/LorentzVector.h"
0005 #include "DataFormats/Math/interface/Vector3D.h"
0006
0007 #include "Math/GenVector/RotationX.h"
0008 #include "Math/GenVector/RotationY.h"
0009 #include "Math/GenVector/RotationZ.h"
0010 #include "Math/GenVector/Rotation3D.h"
0011 #include "Math/GenVector/AxisAngle.h"
0012 #include "Math/GenVector/Boost.h"
0013
0014 #include <string>
0015 #include <iosfwd>
0016
0017 #include <memory>
0018
0019 /**
0020  * A prototype for a particle class.
0021  *  This class describes a general particle beeing a fourvector
0022  *  and containing a vertex (fourvector). It is defined in RawParticle.h
0023  * \author Stephan Wynhoff
0024  */
0025 typedef math::XYZTLorentzVector XYZTLorentzVector;
0026 typedef math::XYZVector XYZVector;
0027
0028 class RawParticle;
0029
0030 namespace rawparticle {
0031   ///Create a particle with momentum 'p' at space-time point xStart
0032   /// The particle will be a muon if iParticle==true, else it will
0033   /// be an anti-muon.
0034   RawParticle makeMuon(bool isParticle, const XYZTLorentzVector& p, const XYZTLorentzVector& xStart);
0035 }  // namespace rawparticle
0036
0037 class RawParticle {
0038 public:
0039   friend RawParticle rawparticle::makeMuon(bool, const XYZTLorentzVector&, const XYZTLorentzVector&);
0040   friend RawParticle unchecked_makeParticle(int id, const math::XYZTLorentzVector& p, double mass, double charge);
0041   friend RawParticle unchecked_makeParticle(
0042       int id, const math::XYZTLorentzVector& p, const math::XYZTLorentzVector& xStart, double mass, double charge);
0043
0044   typedef ROOT::Math::AxisAngle Rotation;
0045   typedef ROOT::Math::Rotation3D Rotation3D;
0046   typedef ROOT::Math::RotationX RotationX;
0047   typedef ROOT::Math::RotationY RotationY;
0048   typedef ROOT::Math::RotationZ RotationZ;
0049   typedef ROOT::Math::Boost Boost;
0050
0051   RawParticle() = default;
0052
0053   /** Construct from a fourvector.
0054    *  The fourvector is taken for the particle, the vertex is set to 0.
0055    */
0056   RawParticle(const XYZTLorentzVector& p);
0057
0058   /** Construct from 2 fourvectors.
0059    *  The first fourvector is taken for the particle, the second for its vertex.
0060    */
0061   RawParticle(const XYZTLorentzVector& p, const XYZTLorentzVector& xStart, double charge = 0.);
0062
0063   /** Construct from fourmomentum components.
0064    *  Vertex is set to 0.
0065    */
0066   RawParticle(double px, double py, double pz, double e, double charge = 0.);
0067
0068   /** Copy constructor    */
0069   RawParticle(const RawParticle& p) = default;
0070   RawParticle(RawParticle&& p) = default;
0071
0072   /** Copy assignment operator */
0073   RawParticle& operator=(const RawParticle& rhs) = default;
0074   RawParticle& operator=(RawParticle&& rhs) = default;
0075
0076 public:
0077   /** Set the status of this particle.
0078    *  The coding follows PYTHIAs convention:
0079    *  1 = stable
0080    */
0081   void setStatus(int istat);
0082
0083   /// set the RECONSTRUCTED mass
0084   void setMass(float m);
0085
0086   /// set the MEASURED charge
0087   void setCharge(float q);
0088
0089   /// set the time of creation
0090   void setT(const double t);
0091
0092   ///  set the vertex
0093   void setVertex(const XYZTLorentzVector& vtx);
0094   void setVertex(double xv, double yv, double zv, double tv);
0095
0096   ///  set the momentum
0097   void setMomentum(const XYZTLorentzVector& vtx);
0098   void setMomentum(double xv, double yv, double zv, double tv);
0099
0100   void SetPx(double);
0101   void SetPy(double);
0102   void SetPz(double);
0103   void SetE(double);
0104
0105   /***  methods to be overloaded to include vertex ***/
0106
0107   /** Boost the particle.
0108    *  The arguments are the \f$\beta\f$ values of the boost in x, y
0109    * and z direction. \warning What happens to the vertex?
0110    */
0111   void boost(double bx, double by, double bz);
0112   void boost(const Boost& b);
0113
0114   //  inline void boost(const Hep3Vector<double> &bv );
0115
0116   /** Rotate the particle around an axis in space.
0117    *  The arguments give the amount to rotate \a rphi in radian and a vector
0118    *  \a raxis in 3D space around which the rotation is done. The vertex is
0119    *  rotated using the same transformation.
0120    */
0121   void rotate(double rphi, const XYZVector& raxis);
0122   void rotate(const Rotation& r);
0123   void rotate(const Rotation3D& r);
0124
0125   /** \warning not yet implemented   */
0126   //   void rotateUz(Hep3Vector &nuz);
0127
0128   /** Rotate around x axis.
0129    *  Rotate \a rphi radian around the x axis. The Vertex is rotated as well.
0130    */
0131   void rotateX(double rphi);
0132   void rotate(const RotationX& r);
0133
0134   /** Rotate around z axis.
0135    *  Rotate \a rphi radian around the z axis. The Vertex is rotated as well.
0136    */
0137
0138   void rotateY(double rphi);
0139   void rotate(const RotationY& r);
0140   /** Rotate around z axis.
0141    *  Rotate \a rphi radian around the z axis. The Vertex is rotated as well.
0142    */
0143
0144   void rotateZ(double rphi);
0145   void rotate(const RotationZ& r);
0146
0147   /** Translate the vertex by a given space amount */
0148   void translate(const XYZVector& t);
0149
0150   //  inline RawParticle & transform(const HepRotation &rot);
0151   //  inline RawParticle & transform(const HepLorentzRotation &rot);
0152
0153   /** Convert the particle to its charge conjugate state.
0154       This operation resets the particle ID to that of the charge conjugated
0155       particle (if one exists). Also the measured charge is multiplied by -1.
0156    */
0157   void chargeConjugate();
0158
0159   int pid() const;  //!< get the HEP particle ID number
0160
0161   int status() const;  //!< get the particle status
0162
0163   double charge() const;  //!< get the MEASURED charge
0164
0165   double mass() const;  //!< get the MEASURED mass
0166
0167   /** Get the pseudo rapidity of the particle.
0168    * \f$ \eta = -\log ( \tan ( \vartheta/2)) \f$
0169    */
0170   double eta() const;
0171
0172   /// Cos**2(theta) is faster to determine than eta
0173   double cos2Theta() const;
0174   double cos2ThetaV() const;
0175
0176   double et() const;  //!< get the transverse energy
0177
0178   double x() const;  //!< x of vertex
0179   double X() const;  //!< x of vertex
0180
0181   double y() const;  //!< y of vertex
0182   double Y() const;  //!< y of vertex
0183
0184   double z() const;  //!< z of vertex
0185   double Z() const;  //!< z of vertex
0186
0187   double t() const;  //!< vertex time
0188   double T() const;  //!< vertex time
0189
0190   double r() const;  //!< vertex radius
0191   double R() const;  //!< vertex radius
0192
0193   double r2() const;  //!< vertex radius**2
0194   double R2() const;  //!< vertex radius**2
0195
0196   const XYZTLorentzVector& vertex() const;  //!< the vertex fourvector
0197
0198   double px() const;  //!< x of the momentum
0199   double Px() const;  //!< x of the momentum
0200
0201   double py() const;  //!< y of the momentum
0202   double Py() const;  //!< y of the momentum
0203
0204   double pz() const;  //!< z of the momentum
0205   double Pz() const;  //!< z of the momentum
0206
0207   double e() const;  //!< energy of the momentum
0208   double E() const;  //!< energy of the momentum
0209
0210   double Pt() const;  //!< transverse momentum
0211   double pt() const;  //!< transverse momentum
0212
0213   double Perp2() const;  //!< perpendicular momentum squared
0214
0215   double mag() const;  //!< the magnitude of the momentum
0216
0217   double theta() const;  //!< theta of momentum vector
0218   double phi() const;    //!< phi of momentum vector
0219
0220   double M2() const;  //!< mass squared
0221
0222   const XYZTLorentzVector& momentum() const;  //!< the momentum fourvector
0223   XYZTLorentzVector& momentum();              //!< the momentum fourvector
0224
0225   XYZVector Vect() const;  //!< the momentum threevector
0226
0227   /** Print the name of the particle.
0228    *  The name is deduced from the particle ID using a particle data table.
0229    *  It is printed with a length of 10 characters. If the id number cannot
0230    *  be found in the table "unknown" is printed as name.
0231    */
0232   void printName() const;
0233
0234   /** Print the formated particle information.
0235    *  The format is:
0236    *  NAME______PX______PY______PZ______E_______Mtheo___Mrec____Qrec____X_______Y_______Z_______T_______
0237    */
0238   void print() const;
0239
0240   /** Is the particle marked as used.
0241    *  The three methods isUsed(), use() and reUse() implement a simple
0242    *  locking mechanism.
0243    */
0244   int isUsed() const { return myUsed; }
0245
0246   /** Lock the particle, see isUsed()
0247    */
0248   void use() { myUsed = 1; }
0249
0250   /** Unlock the particle, see isUsed()
0251    */
0252   void reUse() { myUsed = 0; }
0253
0254 private:
0255   /** Construct from a fourvector and a PID.
0256    *  The fourvector and PID are taken for the particle, the vertex is set to 0.
0257    */
0258   RawParticle(const int id, const XYZTLorentzVector& p, double mass, double charge);
0259
0260   /** Construct from 2 fourvectosr and a PID.
0261    *  The fourvector and PID are taken for the particle, the vertex is set to 0.
0262    */
0263   RawParticle(const int id, const XYZTLorentzVector& p, const XYZTLorentzVector& xStart, double mass, double charge);
0264
0265 private:
0266   XYZTLorentzVector myMomentum;  //!< the four vector of the momentum
0267   XYZTLorentzVector myVertex;    //!< the four vector of the vertex
0268   double myCharge = 0.;          //!< the MEASURED charge
0269   double myMass = 0.;            //!< the RECONSTRUCTED mass
0270   int myId = 0;                  //!< the particle id number HEP-PID
0271   int myStatus = 99;             //!< the status code according to PYTHIA
0272   int myUsed = 0;                //!< status of the locking
0273 };
0274
0275 std::ostream& operator<<(std::ostream& o, const RawParticle& p);
0276
0277 inline int RawParticle::pid() const { return myId; }
0278 inline int RawParticle::status() const { return myStatus; }
0279 inline double RawParticle::eta() const { return -std::log(std::tan(this->theta() / 2.)); }
0280 inline double RawParticle::cos2Theta() const { return Pz() * Pz() / myMomentum.Vect().Mag2(); }
0281 inline double RawParticle::cos2ThetaV() const { return Z() * Z() / myVertex.Vect().Mag2(); }
0282 inline double RawParticle::x() const { return myVertex.Px(); }
0283 inline double RawParticle::y() const { return myVertex.Py(); }
0284 inline double RawParticle::z() const { return myVertex.Pz(); }
0285 inline double RawParticle::t() const { return myVertex.E(); }
0286 inline double RawParticle::X() const { return myVertex.Px(); }
0287 inline double RawParticle::Y() const { return myVertex.Py(); }
0288 inline double RawParticle::Z() const { return myVertex.Pz(); }
0289 inline double RawParticle::T() const { return myVertex.E(); }
0290 inline double RawParticle::R() const { return std::sqrt(R2()); }
0291 inline double RawParticle::R2() const { return myVertex.Perp2(); }
0292 inline double RawParticle::r() const { return std::sqrt(r2()); }
0293 inline double RawParticle::r2() const { return myVertex.Perp2(); }
0294 inline double RawParticle::charge() const { return myCharge; }
0295 inline double RawParticle::mass() const { return myMass; }
0296 inline double RawParticle::px() const { return myMomentum.px(); }
0297 inline double RawParticle::Px() const { return myMomentum.Px(); }
0298
0299 inline double RawParticle::py() const { return myMomentum.py(); }
0300 inline double RawParticle::Py() const { return myMomentum.Py(); }
0301
0302 inline double RawParticle::pz() const { return myMomentum.pz(); }
0303 inline double RawParticle::Pz() const { return myMomentum.Pz(); }
0304
0305 inline double RawParticle::e() const { return myMomentum.e(); }
0306 inline double RawParticle::E() const { return myMomentum.E(); }
0307
0308 inline double RawParticle::Pt() const { return myMomentum.Pt(); }
0309 inline double RawParticle::pt() const { return myMomentum.pt(); }
0310
0311 inline double RawParticle::Perp2() const { return myMomentum.Perp2(); }
0312
0313 inline double RawParticle::mag() const { return myMomentum.mag(); }
0314
0315 inline double RawParticle::theta() const { return myMomentum.theta(); }
0316 inline double RawParticle::phi() const { return myMomentum.phi(); }
0317
0318 inline double RawParticle::M2() const { return myMomentum.M2(); }
0319
0320 inline const XYZTLorentzVector& RawParticle::vertex() const { return myVertex; }
0321 inline const XYZTLorentzVector& RawParticle::momentum() const { return myMomentum; }
0322 inline XYZTLorentzVector& RawParticle::momentum() { return myMomentum; }
0323 inline XYZVector RawParticle::Vect() const { return myMomentum.Vect(); }
0324
0325 inline void RawParticle::setVertex(const XYZTLorentzVector& vtx) { myVertex = vtx; }
0326 inline void RawParticle::setVertex(double a, double b, double c, double d) { myVertex.SetXYZT(a, b, c, d); }
0327
0328 inline void RawParticle::setMomentum(const XYZTLorentzVector& p4) { myMomentum = p4; }
0329 inline void RawParticle::setMomentum(double a, double b, double c, double d) { myMomentum.SetXYZT(a, b, c, d); }
0330
0331 inline void RawParticle::SetPx(double px) { myMomentum.SetPx(px); }
0332 inline void RawParticle::SetPy(double py) { myMomentum.SetPy(py); }
0333 inline void RawParticle::SetPz(double pz) { myMomentum.SetPz(pz); }
0334 inline void RawParticle::SetE(double e) { myMomentum.SetE(e); }
0335
0336 inline void RawParticle::rotate(const RawParticle::Rotation3D& r) {
0337   XYZVector v(r(myMomentum.Vect()));
0338   setMomentum(v.X(), v.Y(), v.Z(), E());
0339 }
0340
0341 inline void RawParticle::rotate(const RawParticle::Rotation& r) {
0342   XYZVector v(r(myMomentum.Vect()));
0343   setMomentum(v.X(), v.Y(), v.Z(), E());
0344 }
0345
0346 inline void RawParticle::rotate(const RawParticle::RotationX& r) {
0347   XYZVector v(r(myMomentum.Vect()));
0348   setMomentum(v.X(), v.Y(), v.Z(), E());
0349 }
0350
0351 inline void RawParticle::rotate(const RawParticle::RotationY& r) {
0352   XYZVector v(r(myMomentum.Vect()));
0353   setMomentum(v.X(), v.Y(), v.Z(), E());
0354 }
0355
0356 inline void RawParticle::rotate(const RawParticle::RotationZ& r) {
0357   XYZVector v(r(myMomentum.Vect()));
0358   setMomentum(v.X(), v.Y(), v.Z(), E());
0359 }
0360
0361 inline void RawParticle::boost(const RawParticle::Boost& b) {
0362   XYZTLorentzVector p(b(momentum()));
0363   setMomentum(p.Px(), p.Py(), p.Pz(), p.E());
0364 }
0365
0366 inline void RawParticle::translate(const XYZVector& tr) {
0367   myVertex.SetXYZT(X() + tr.X(), Y() + tr.Y(), Z() + tr.Z(), T());
0368 }
0369
0370 #endif