src/framework/FourVector.h

0001 /*******************************************************************************
0002  * Copyright (c) The JETSCAPE Collaboration, 2018
0003  *
0004  * Modular, task-based framework for simulating all aspects of heavy-ion collisions
0005  *
0006  * For the list of contributors see AUTHORS.
0007  *
0008  * Report issues at https://github.com/JETSCAPE/JETSCAPE/issues
0009  *
0010  * or via email to bugs.jetscape@gmail.com
0011  *
0012  * Distributed under the GNU General Public License 3.0 (GPLv3 or later).
0013  * See COPYING for details.
0014  ******************************************************************************/
0015
0016 #ifndef FOURVECTOR_H
0017 #define FOURVECTOR_H
0018
0019 #include <iostream>
0020 #include <cmath>
0021 #include <vector>
0022 #include <cstdlib>
0023 #include <climits>
0024
0025 using std::cout;
0026 using std::cerr;
0027 using std::endl;
0028
0029 namespace Jetscape {
0030
0031 class FourVector {
0032   // the class of four vectors
0033
0034 public:
0035   FourVector() //default constructor
0036   {
0037     tv = xv = yv = zv = 0.0;
0038   };
0039
0040   FourVector(const FourVector &srv)
0041       : xv(srv.xv), yv(srv.yv), zv(srv.zv), tv(srv.tv){}; // copy constructor
0042
0043   FourVector(double a[4]) // constructor with array input
0044   {
0045     tv = a[0];
0046     xv = a[1];
0047     yv = a[2];
0048     zv = a[3];
0049   };
0050
0051   FourVector(double x_in, double y_in, double z_in,
0052              double t_in) // constructor with component input
0053   {
0054     tv = t_in;
0055     xv = x_in;
0056     yv = y_in;
0057     zv = z_in;
0058   };
0059
0060   void clear() { tv = xv = yv = zv = 0.0; }
0061
0062   // constructors do all sets
0063
0064   void Set(double x_in, double y_in, double z_in, double t_in) {
0065     tv = t_in;
0066     xv = x_in;
0067     yv = y_in;
0068     zv = z_in;
0069   }
0070
0071   void Set(double a[4]) {
0072     tv = a[0];
0073     xv = a[1];
0074     yv = a[2];
0075     zv = a[3];
0076   };
0077
0078   // all gets are done with name calls e.g., vec.x()
0079   double x() const { return (xv); };
0080
0081   double y() const { return (yv); };
0082
0083   double z() const { return (zv); };
0084
0085   double t() const { return (tv); };
0086
0087   const double comp(int i) const {
0088     switch (i) {
0089     case 0:
0090       return (tv);
0091       break;
0092     case 1:
0093       return (xv);
0094       break;
0095     case 2:
0096       return (yv);
0097       break;
0098     case 3:
0099       return (zv);
0100       break;
0101     default:
0102       cout << " component index beyond 0-3! Returning garbage ..." << endl;
0103       return (a_very_large_number);
0104       break;
0105     }
0106   }
0107
0108   double plus() { return ((zv + tv) / sqrt(2.0)); };
0109
0110   double minus() { return ((tv - zv) / sqrt(2.0)); };
0111
0112   double rapidity() {
0113     if (this->minus() > 0.0)
0114       return (std::log(this->plus() / this->minus()) / 2.0);
0115     cout << endl
0116          << "ERROR: z component exceeds t component, cannot calculate rapidity"
0117          << endl;
0118     return (0);
0119   };
0120
0121   double phi() {
0122     if (fabs(x()) < rounding_error && fabs(y()) < rounding_error) {
0123       return 0;
0124     }
0125     double phi = atan2(y(), x());
0126     while (phi < 0)
0127       phi += 2.0 * pi;
0128     return phi;
0129   };
0130
0131   double operator*(FourVector &c) {
0132     return (tv * c.t() - xv * c.x() - yv * c.y() - zv * c.z());
0133   };
0134
0135   FourVector &operator+=(FourVector &c) {
0136     tv += c.t();
0137     xv += c.x();
0138     yv += c.y();
0139     zv += c.z();
0140
0141     return (*this);
0142   };
0143
0144   FourVector &operator-=(FourVector &c) {
0145     tv -= c.t();
0146     xv -= c.x();
0147     yv -= c.y();
0148     zv -= c.z();
0149
0150     return (*this);
0151   };
0152
0153   FourVector &operator=(FourVector &c) {
0154     tv = c.t();
0155     xv = c.x();
0156     yv = c.y();
0157     zv = c.z();
0158     return (*this);
0159   };
0160
0161   FourVector &operator=(const FourVector &c) {
0162     tv = c.tv;
0163     xv = c.xv;
0164     yv = c.yv;
0165     zv = c.zv;
0166     return (*this);
0167   };
0168
0169   void rotate_around_z(double theta) {
0170     double new_xv, new_yv;
0171
0172     new_xv = xv * cos(theta) - yv * sin(theta);
0173
0174     new_yv = yv * cos(theta) + xv * sin(theta);
0175
0176     xv = new_xv;
0177     yv = new_yv;
0178   };
0179
0180 private:
0181   // the v is for vector, we call the private variables, xv, tv etc., so that get function
0182   // calls will be called x, t etc.
0183   double xv, yv, zv, tv;
0184 };
0185
0186 }; // namespace Jetscape
0187
0188 #endif // FOURVECTOR_H