Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:49:54

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