Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
rmn.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file rmn.cpp
1 #include <iostream>
2 #include <fstream>
3 #include <iomanip>
4 #include <cstdlib>
5 #include <execinfo.h>
6 #include <signal.h>
7 #include "rmn.h"
8 
9 using namespace std;
10 
11 // flag optionally needed to track
12 // problems in the transformation procedures
13 bool debugRiemann;
14 
15 void handler(int sig) {
16  void *array[10];
17  size_t size;
18 
19  // get void*'s for all entries on the stack
20  size = backtrace(array, 10);
21 
22  // print out all the frames to stderr
23  cout << "Error: signal " << sig << endl;
24  backtrace_symbols_fd(array, size, 2);
25  exit(1);
26 }
27 
28 void transformPV(EoS *eos, double Q[7], double &e, double &p, double &nb,
29  double &nq, double &ns, double &vx, double &vy, double &vz) {
30  // conserved -> primitive transtormation requires
31  // a numerical solution to 1D nonlinear algebraic equation:
32  // v = M / ( Q_t + p(Q_t-M*v, n) ) (A.2)
33  // M being the modulo of the vector {Q_x, Q_y, Q_z}.
34  // Bisection/Newton methods are used to solve the equation.
35  const int MAXIT = 100; // maximum number of iterations
36  const double dpe = 1. / 3.; // dp/de estimate for Newton method
37  const double corrf = 0.9999; // corrected value of M
38  // when it brakes the speed of light limit, M>Q_t
39  double v, vl = 0., vh = 1., dvold, dv, f, df;
40  if (debugRiemann) {
41  cout << "transformPV debug---------------\n";
42  cout << setw(14) << Q[0] << setw(14) << Q[1] << setw(14) << Q[2] << setw(14)
43  << Q[3] << endl;
44  cout << setw(14) << Q[4] << setw(14) << Q[5] << setw(14) << Q[6] << endl;
45  }
46  double M = sqrt(Q[X_] * Q[X_] + Q[Y_] * Q[Y_] + Q[Z_] * Q[Z_]);
47  if (Q[T_] <= 0.) {
48  e = 0.;
49  p = 0.;
50  vx = vy = vz = 0.;
51  nb = nq = ns = 0.;
52  return;
53  }
54  if (M == 0.) {
55  e = Q[T_];
56  vx = 0.;
57  vy = 0.;
58  vz = 0.;
59  nb = Q[NB_];
60  nq = Q[NQ_];
61  ns = Q[NS_];
62  p = eos->p(e, nb, nq, ns);
63  return;
64  }
65  if (M > Q[T_]) {
66  Q[X_] *= corrf * Q[T_] / M;
67  Q[Y_] *= corrf * Q[T_] / M;
68  Q[Z_] *= corrf * Q[T_] / M;
69  M = Q[T_] * corrf;
70  }
71 
72  v = 0.5 * (vl + vh);
73  e = Q[T_] - M * v;
74  if (e < 0.) e = 0.;
75  nb = Q[NB_] * sqrt(1 - v * v);
76  nq = Q[NQ_] * sqrt(1 - v * v);
77  ns = Q[NS_] * sqrt(1 - v * v);
78  p = eos->p(e, nb, nq, ns);
79  f = (Q[T_] + p) * v - M;
80  df = (Q[T_] + p) - M * v * dpe;
81  dvold = vh - vl;
82  dv = dvold;
83  for (int i = 0; i < MAXIT; i++) {
84  if ((f + df * (vh - v)) * (f + df * (vl - v)) > 0. ||
85  fabs(2. * f) > fabs(dvold * df)) { // bisection
86  dvold = dv;
87  dv = 0.5 * (vh - vl);
88  v = vl + dv;
89  // cout << "BISECTION v = " << setw(12) << v << endl ;
90  } else { // Newton
91  dvold = dv;
92  dv = f / df;
93  v -= dv;
94  // cout << "NEWTON v = " << setw(12) << v << endl ;
95  }
96  if (fabs(dv) < 0.00001) break;
97 
98  e = Q[T_] - M * v;
99  if (e < 0.) e = 0.;
100  nb = Q[NB_] * sqrt(1 - v * v);
101  nq = Q[NQ_] * sqrt(1 - v * v);
102  ns = Q[NS_] * sqrt(1 - v * v);
103  p = eos->p(e, nb, nq, ns);
104  f = (Q[T_] + p) * v - M;
105  df = (Q[T_] + p) - M * v * dpe;
106 
107  if (f > 0.)
108  vh = v;
109  else
110  vl = v;
111  if (debugRiemann)
112  cout << "step " << i << " " << e << " " << nb << " " << nq << " "
113  << ns << " " << p << " " << vx << " " << vy << " " << vz << endl;
114  // if(i==40) { cout << "error : does not converge\n" ; exit(1) ; } ;
115  } // for loop
116  //----------after
117  // v = 0.5*(vh+vl) ;
118  vx = v * Q[X_] / M;
119  vy = v * Q[Y_] / M;
120  vz = v * Q[Z_] / M;
121  e = Q[T_] - M * v;
122  p = eos->p(e, nb, nq, ns);
123  nb = Q[NB_] * sqrt(1 - vx * vx - vy * vy - vz * vz);
124  nq = Q[NQ_] * sqrt(1 - vx * vx - vy * vy - vz * vz);
125  ns = Q[NS_] * sqrt(1 - vx * vx - vy * vy - vz * vz);
126  if (e < 0. || sqrt(vx * vx + vy * vy + vz * vz) > 1.) {
127  cout << Q[T_] << " " << Q[X_] << " " << Q[Y_] << " " << Q[Z_] << " "
128  << Q[NB_] << endl;
129  cout << "transformRF::Error\n";
130  }
131  if (!(nb < 0. || nb >= 0.)) {
132  cout << "transformRF::Error nb=#ind\n";
133  // return ;
134  }
135 }
136 
137 void transformPVBulk(EoS *eos, double Pi, double Q[7], double &e, double &p,
138  double &nb, double &nq, double &ns, double &vx, double &vy,
139  double &vz) {
140  const int MAXIT = 100;
141  const double dpe = 1. / 3.;
142  const double corrf = 0.9999;
143  double v, vl = 0., vh = 1., dvold, dv, f, df;
144  if (debugRiemann) {
145  cout << "transformPVBulk debug---------------\n";
146  cout << setw(14) << Q[0] << setw(14) << Q[1] << setw(14) << Q[2] << setw(14)
147  << Q[3] << endl;
148  cout << setw(14) << Q[4] << setw(14) << Q[5] << setw(14) << Q[6] << endl;
149  }
150  double M = sqrt(Q[X_] * Q[X_] + Q[Y_] * Q[Y_] + Q[Z_] * Q[Z_]);
151  if (Q[T_] <= 0.) {
152  e = 0.;
153  p = 0.;
154  vx = vy = vz = 0.;
155  nb = nq = ns = 0.;
156  return;
157  }
158  if (M == 0.) {
159  e = Q[T_];
160  vx = 0.;
161  vy = 0.;
162  vz = 0.;
163  nb = Q[NB_];
164  nq = Q[NQ_];
165  ns = Q[NS_];
166  p = eos->p(e, nb, nq, ns) + Pi;
167  return;
168  }
169  if (M > Q[T_]) {
170  Q[X_] *= corrf * Q[T_] / M;
171  Q[Y_] *= corrf * Q[T_] / M;
172  Q[Z_] *= corrf * Q[T_] / M;
173  M = Q[T_] * corrf;
174  }
175 
176  v = 0.5 * (vl + vh);
177  e = Q[T_] - M * v;
178  if (e < 0.) e = 0.;
179  nb = Q[NB_] * sqrt(1 - v * v);
180  nq = Q[NQ_] * sqrt(1 - v * v);
181  ns = Q[NS_] * sqrt(1 - v * v);
182  p = eos->p(e, nb, nq, ns) + Pi;
183  f = (Q[T_] + p) * v - M;
184  df = (Q[T_] + p) - M * v * dpe;
185  dvold = vh - vl;
186  dv = dvold;
187  for (int i = 0; i < MAXIT; i++) {
188  if ((f + df * (vh - v)) * (f + df * (vl - v)) > 0. ||
189  fabs(2. * f) > fabs(dvold * df)) { // bisection
190  dvold = dv;
191  dv = 0.5 * (vh - vl);
192  v = vl + dv;
193  // cout << "BISECTION v = " << setw(12) << v << endl ;
194  } else { // Newton
195  dvold = dv;
196  dv = f / df;
197  v -= dv;
198  // cout << "NEWTON v = " << setw(12) << v << endl ;
199  }
200  if (fabs(dv) < 0.00001) break;
201 
202  e = Q[T_] - M * v;
203  if (e < 0.) e = 0.;
204  nb = Q[NB_] * sqrt(1 - v * v);
205  nq = Q[NQ_] * sqrt(1 - v * v);
206  ns = Q[NS_] * sqrt(1 - v * v);
207  p = eos->p(e, nb, nq, ns) + Pi;
208  f = (Q[T_] + p) * v - M;
209  df = (Q[T_] + p) - M * v * dpe;
210 
211  if (f > 0.)
212  vh = v;
213  else
214  vl = v;
215  if (debugRiemann)
216  cout << "step " << i << " " << e << " " << nb << " " << nq << " "
217  << ns << " " << p << " " << vx << " " << vy << " " << vz << endl;
218  // if(i==40) { cout << "error : does not converge\n" ; exit(1) ; } ;
219  } // for loop
220  //----------after
221  // v = 0.5*(vh+vl) ;
222  vx = v * Q[X_] / M;
223  vy = v * Q[Y_] / M;
224  vz = v * Q[Z_] / M;
225  e = Q[T_] - M * v;
226  p = eos->p(e, nb, nq, ns);
227  nb = Q[NB_] * sqrt(1 - vx * vx - vy * vy - vz * vz);
228  nq = Q[NQ_] * sqrt(1 - vx * vx - vy * vy - vz * vz);
229  ns = Q[NS_] * sqrt(1 - vx * vx - vy * vy - vz * vz);
230  if (e < 0. || sqrt(vx * vx + vy * vy + vz * vz) > 1.) {
231  cout << Q[T_] << " " << Q[X_] << " " << Q[Y_] << " " << Q[Z_] << " "
232  << Q[NB_] << endl;
233  cout << "transformRF::Error\n";
234  }
235  if (!(nb < 0. || nb >= 0.)) {
236  cout << "transformRF::Error nb=#ind\n";
237  cout << "Q [7]: " << Q[0] << " " << Q[1] << " " << Q[2] << " " << Q[3]
238  << " " << Q[4] << " " << Q[5] << " " << Q[6] << endl;
239  cout << "e vx vy vz nb nq ns: " << e << " " << vx << " " << vy << " " << vz
240  << " " << nb << " " << nq << " " << ns << endl;
241  handler(333);
242  // return ;
243  }
244 }
245 
246 void transformCV(double e, double p, double nb, double nq, double ns, double vx,
247  double vy, double vz, double Q[]) {
248  const double gamma2 = 1. / (1 - vx * vx - vy * vy - vz * vz);
249  Q[T_] = (e + p * (vx * vx + vy * vy + vz * vz)) * gamma2;
250  Q[X_] = (e + p) * vx * gamma2;
251  Q[Y_] = (e + p) * vy * gamma2;
252  Q[Z_] = (e + p) * vz * gamma2;
253  Q[NB_] = nb * sqrt(gamma2);
254  Q[NQ_] = nq * sqrt(gamma2);
255  Q[NS_] = ns * sqrt(gamma2);
256 }