Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4Utils.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4Utils.cc
1 #include "PHG4Utils.h"
2 
3 #include <phool/phool.h>
4 
5 #include <Geant4/G4Colour.hh> // for G4Colour
6 #include <Geant4/G4VisAttributes.hh>
7 
8 #include <TSystem.h>
9 
10 #pragma GCC diagnostic push
11 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
12 #include <boost/algorithm/string.hpp>
13 #pragma GCC diagnostic pop
14 
15 #include <boost/algorithm/hex.hpp>
16 #include <boost/uuid/detail/md5.hpp>
17 
18 #include <algorithm> // for copy
19 #include <cmath>
20 #include <cstdlib> // for exit
21 #include <fstream>
22 #include <iostream> // for operator<<, endl, basic_ostream
23 #include <iterator> // for back_insert_iterator
24 #include <vector> // for vector
25 
26 using namespace std;
27 
28 double PHG4Utils::_eta_coverage = 1.;
29 
30 double
31 PHG4Utils::GetLengthForRapidityCoverage(const double radius, const double eta)
32 {
33  double length;
34  double theta = 2.0 * std::atan(std::exp(-eta));
35  length = radius / std::tan(theta);
36  return length;
37 }
38 
39 double
41 {
42  return GetLengthForRapidityCoverage(radius, _eta_coverage);
43 }
44 
46 {
47  _eta_coverage = eta;
48 }
49 
50 double
52 {
53  double theta = 2 * atan(exp(-eta));
54  return theta;
55 }
56 
57 double
59 {
60  double eta = -log(tan(theta / 2.));
61  return eta;
62 }
63 
64 pair<double, double>
65 PHG4Utils::get_etaphi(const double x, const double y, const double z)
66 {
67  double eta;
68  double phi;
69  double radius;
70  double theta;
71  radius = sqrt(x * x + y * y);
72  phi = atan2(y, x);
73  theta = atan2(radius, z);
74  eta = -log(tan(theta / 2.));
75  return make_pair(eta, phi);
76 }
77 
78 double
79 PHG4Utils::get_eta(const double radius, const double z)
80 {
81  double eta;
82  double theta;
83  theta = atan2(radius, fabs(z));
84  eta = -log(tan(theta / 2.));
85  if (z < 0)
86  {
87  eta = -eta;
88  }
89  return eta;
90 }
91 
92 void PHG4Utils::SetColour(G4VisAttributes* att, const string& material)
93 {
94  if (!att)
95  {
96  cout << "G4VisAttributes pointer is NULL" << endl;
97  return;
98  }
99  if (material == "AL_BABAR_MAG")
100  {
101  att->SetColour(G4Colour::Blue());
102  }
103  else if (material == "BlackHole")
104  {
105  att->SetColour(G4Colour::Black());
106  }
107  else if (material == "C4F10")
108  {
109  att->SetColour(0., 0., 0.5, 0.25);
110  }
111  else if (material == "CF4")
112  {
113  att->SetColour(G4Colour::Magenta());
114  }
115  else if (material == "G4_AIR")
116  {
117  att->SetColour(G4Colour::Black());
118  }
119  else if (material == "G4_Al")
120  {
121  att->SetColour(G4Colour::Gray());
122  }
123  else if (material == "G4_Au")
124  {
125  att->SetColour(G4Colour::Yellow());
126  }
127  else if (material == "G4_CARBON_DIOXIDE")
128  {
129  att->SetColour(G4Colour::Green());
130  }
131  else if (material == "G4_CELLULOSE_CELLOPHANE")
132  {
133  att->SetColour(0.25, 0.25, 0.);
134  }
135  else if (material == "G4_Cu")
136  {
137  att->SetColour(1., 0.51, 0.278);
138  }
139  else if (material == "G4_Fe")
140  {
141  att->SetColour(0.29, 0.44, 0.54);
142  }
143  else if (material == "G4_KAPTON")
144  {
145  att->SetColour(G4Colour::Yellow());
146  }
147  else if (material == "G4_MYLAR")
148  {
149  att->SetColour(0.5, 0.5, 0.5, 0.25);
150  }
151  else if (material == "G4_METHANE")
152  {
153  att->SetColour(0., 1., 1., 0.25);
154  }
155  else if (material == "G4_Si")
156  {
157  att->SetColour(G4Colour::Yellow());
158  }
159  else if (material == "G4_TEFLON")
160  {
161  att->SetColour(G4Colour::White());
162  }
163  else if (material == "G4_W")
164  {
165  att->SetColour(0.36, 0.36, 0.36);
166  }
167  else if (material == "Quartz")
168  {
169  att->SetColour(G4Colour::Green());
170  }
171  else if (material == "Scintillator" || material == "G4_POLYSTYRENE")
172  {
173  att->SetColour(0., 1., 1.);
174  }
175  else if (material == "W_Epoxy")
176  {
177  att->SetColour(0.5, 0.5, 0.5);
178  }
179  else if (material == "G10")
180  {
181  att->SetColour(1., 1., 0., 0.5);
182  }
183  else
184  {
185  //cout << "default color red for material " << material << endl;
186  att->SetColour(G4Colour::Cyan());
187  }
188  return;
189 }
190 // returns success/failure and intersection point (x/y)
191 std::pair<bool, std::pair<double, double>> PHG4Utils::lines_intersect(
192  double ax,
193  double ay,
194  double bx,
195  double by,
196  double cx,
197  double cy,
198  double dx,
199  double dy)
200 {
201  // Find if a line segment limited by points A and B
202  // intersects line segment limited by points C and D.
203  // First check if an infinite line defined by A and B intersects
204  // segment (C,D). If h is from 0 to 1 line and line segment intersect
205  // Then check in intersection point is between C and D
206  double ex = bx - ax; // E=B-A
207  double ey = by - ay;
208  double fx = dx - cx; // F=D-C
209  double fy = dy - cy;
210  double px = -ey; // P
211  double py = ex;
212 
213  double bottom = fx * px + fy * py; // F*P
214  double gx = ax - cx; // A-C
215  double gy = ay - cy;
216  double top = gx * px + gy * py; // G*P
217 
218  double h = 99999.;
219  if (bottom != 0.)
220  {
221  h = top / bottom;
222  }
223 
224  //intersection point R = C + F*h
225  if (h > 0. && h < 1.)
226  {
227  double rx = cx + fx * h;
228  double ry = cy + fy * h;
229  //cout << " line/segment intersection coordinates: " << *rx << " " << *ry << endl;
230  if ((rx > ax && rx > bx) || (rx < ax && rx < bx) || (ry < ay && ry < by) || (ry > ay && ry > by))
231  {
232  //cout << " NO segment/segment intersection!" << endl;
233  return make_pair(false, make_pair(NAN, NAN));
234  }
235  else
236  {
237  //cout << " segment/segment intersection!" << endl;
238  return make_pair(true, make_pair(rx, ry));
239  }
240  }
241 
242  return make_pair(false, make_pair(NAN, NAN));
243 }
244 
245 // returns success/failure and length of the line segment inside the rectangle (output)
247  double ax,
248  double ay,
249  double bx,
250  double by,
251  double cx,
252  double cy,
253  double dx,
254  double dy)
255 {
256  // find if a line isegment limited by points (A,B)
257  // intersects with a rectangle defined by two
258  // corner points (C,D) two other points are E and F
259  // E--------D
260  // | |
261  // | |
262  // C--------F
263 
264  if (cx > dx || cy > dy)
265  {
266  cout << PHWHERE << "ERROR: Bad rectangle definition!" << endl;
267  return make_pair(false, NAN);
268  }
269 
270  double ex = cx;
271  double ey = dy;
272  double fx = dx;
273  double fy = cy;
274 
275  vector<double> vx;
276  vector<double> vy;
277  pair<bool, pair<double, double>> intersect1 = lines_intersect(ax, ay, bx, by, cx, cy, fx, fy);
278  // bool i1 = lines_intersect(ax, ay, bx, by, cx, cy, fx, fy, &rx, &ry);
279  if (intersect1.first)
280  {
281  vx.push_back(intersect1.second.first);
282  vy.push_back(intersect1.second.second);
283  }
284  pair<bool, pair<double, double>> intersect2 = lines_intersect(ax, ay, bx, by, fx, fy, dx, dy);
285  // bool i2 = lines_intersect(ax, ay, bx, by, fx, fy, dx, dy, &rx, &ry);
286  if (intersect2.first)
287  {
288  vx.push_back(intersect2.second.first);
289  vy.push_back(intersect2.second.second);
290  }
291  pair<bool, pair<double, double>> intersect3 = lines_intersect(ax, ay, bx, by, ex, ey, dx, dy);
292  // bool i3 = lines_intersect(ax, ay, bx, by, ex, ey, dx, dy, &rx, &ry);
293  if (intersect3.first)
294  {
295  vx.push_back(intersect3.second.first);
296  vy.push_back(intersect3.second.second);
297  }
298  pair<bool, pair<double, double>> intersect4 = lines_intersect(ax, ay, bx, by, cx, cy, ex, ey);
299  // bool i4 = lines_intersect(ax, ay, bx, by, cx, cy, ex, ey, &rx, &ry);
300  if (intersect4.first)
301  {
302  vx.push_back(intersect4.second.first);
303  vy.push_back(intersect4.second.second);
304  }
305 
306  //cout << "Rectangle intersections: " << i1 << " " << i2 << " " << i3 << " " << i4 << endl;
307  //cout << "Number of intersections = " << vx.size() << endl;
308 
309  double rr = 0.;
310  if (vx.size() == 2)
311  {
312  rr = sqrt((vx[0] - vx[1]) * (vx[0] - vx[1]) + (vy[0] - vy[1]) * (vy[0] - vy[1]));
313  // cout << "Length of intersection = " << *rr << endl;
314  }
315  if (vx.size() == 1)
316  {
317  // find which point (A or B) is within the rectangle
318  if (ax > cx && ay > cy && ax < dx && ay < dy) // point A is inside the rectangle
319  {
320  //cout << "Point A is inside the rectangle." << endl;
321  rr = sqrt((vx[0] - ax) * (vx[0] - ax) + (vy[0] - ay) * (vy[0] - ay));
322  }
323  if (bx > cx && by > cy && bx < dx && by < dy) // point B is inside the rectangle
324  {
325  //cout << "Point B is inside the rectangle." << endl;
326  rr = sqrt((vx[0] - bx) * (vx[0] - bx) + (vy[0] - by) * (vy[0] - by));
327  }
328  }
329 
330  if (intersect1.first || intersect2.first || intersect3.first || intersect4.first)
331  {
332  return make_pair(true, rr);
333  }
334  return make_pair(false, NAN);
335 }
336 
337 double PHG4Utils::sA(double r, double x, double y)
338 {
339  // Uses analytic formula for the integral of a circle between limits set by the corner of a rectangle
340  // It is called repeatedly to find the overlap area between the circle and rectangle
341  // I found this code implementing the integral on a web forum called "ars technica",
342  // https://arstechnica.com/civis/viewtopic.php?t=306492
343  // posted by "memp"
344 
345  double a;
346 
347  if (x < 0)
348  {
349  return -sA(r, -x, y);
350  }
351 
352  if (y < 0)
353  {
354  return -sA(r, x, -y);
355  }
356 
357  if (x > r)
358  {
359  x = r;
360  }
361 
362  if (y > r)
363  {
364  y = r;
365  }
366 
367  if (x * x + y * y > r * r)
368  {
369  a = r * r * asin(x / r) + x * sqrt(r * r - x * x) + r * r * asin(y / r) + y * sqrt(r * r - y * y) - r * r * M_PI_2;
370 
371  a *= 0.5;
372  }
373  else
374  {
375  a = x * y;
376  }
377 
378  return a;
379 }
380 
381 double PHG4Utils::circle_rectangle_intersection(double x1, double y1, double x2, double y2, double mx, double my, double r)
382 {
383  // Find the area of overlap of a circle and rectangle
384  // Calls sA, which uses an analytic formula to determine the integral of the circle between limits set by the corners of the rectangle
385 
386  // move the rectangle to the frame where the circle is at (0,0)
387  x1 -= mx;
388  x2 -= mx;
389  y1 -= my;
390  y2 -= my;
391 
392  // {
393  // cout << " mx " << mx << " my " << my << " r " << r << " x1 " << x1 << " x2 " << x2 << " y1 " << y1 << " y2 " << y2 << endl;
394  // cout << " sA21 " << sA(r, x2, y1)
395  // << " sA11 " << sA(r, x1, y1)
396  // << " sA22 " << sA(r, x2, y2)
397  // << " sA12 " << sA(r, x1, y2)
398  // << endl;
399  // }
400 
401  return sA(r, x2, y1) - sA(r, x1, y1) - sA(r, x2, y2) + sA(r, x1, y2);
402 }
403 
405 {
406  ifstream myfile;
407  myfile.open(filename);
408  if (!myfile.is_open())
409  {
410  cout << "Error opening " << filename << endl;
411  gSystem->Exit(1);
412  exit(1);
413  }
414  boost::uuids::detail::md5 hash;
415  boost::uuids::detail::md5::digest_type digest;
416  char c;
417  while (myfile.get(c))
418  {
419  hash.process_bytes(&c, 1);
420  }
421  myfile.close();
422  hash.get_digest(digest);
423  const auto charDigest = reinterpret_cast<const char*>(&digest);
424  std::string result;
425  boost::algorithm::hex(charDigest, charDigest + sizeof(boost::uuids::detail::md5::digest_type), std::back_inserter(result));
426  boost::algorithm::to_lower(result);
427  return result;
428 }