g4simulation/g4main/PHG4Utils.cc

0001 #include "PHG4Utils.h"
0002
0003 #include <phool/phool.h>
0004
0005 #include <Geant4/G4Colour.hh>  // for G4Colour
0006 #include <Geant4/G4VisAttributes.hh>
0007
0008 #include <TSystem.h>
0009
0010 #pragma GCC diagnostic push
0011 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0012 #include <boost/algorithm/string.hpp>
0013 #pragma GCC diagnostic pop
0014
0015 #include <boost/algorithm/hex.hpp>
0016 #include <boost/uuid/detail/md5.hpp>
0017
0018 #include <algorithm>  // for copy
0019 #include <cmath>
0020 #include <cstdlib>  // for exit
0021 #include <fstream>
0022 #include <iostream>  // for operator<<, endl, basic_ostream
0023 #include <iterator>  // for back_insert_iterator
0024 #include <vector>    // for vector
0025
0026 using namespace std;
0027
0028 double PHG4Utils::_eta_coverage = 1.;
0029
0030 double
0031 PHG4Utils::GetLengthForRapidityCoverage(const double radius, const double eta)
0032 {
0033   double length;
0034   double theta = 2.0 * std::atan(std::exp(-eta));
0035   length = radius / std::tan(theta);
0036   return length;
0037 }
0038
0039 double
0040 PHG4Utils::GetLengthForRapidityCoverage(const double radius)
0041 {
0042   return GetLengthForRapidityCoverage(radius, _eta_coverage);
0043 }
0044
0045 void PHG4Utils::SetPseudoRapidityCoverage(const double eta)
0046 {
0047   _eta_coverage = eta;
0048 }
0049
0050 double
0051 PHG4Utils::get_theta(const double eta)
0052 {
0053   double theta = 2 * atan(exp(-eta));
0054   return theta;
0055 }
0056
0057 double
0058 PHG4Utils::get_eta(const double theta)
0059 {
0060   double eta = -log(tan(theta / 2.));
0061   return eta;
0062 }
0063
0064 pair<double, double>
0065 PHG4Utils::get_etaphi(const double x, const double y, const double z)
0066 {
0067   double eta;
0068   double phi;
0069   double radius;
0070   double theta;
0071   radius = sqrt(x * x + y * y);
0072   phi = atan2(y, x);
0073   theta = atan2(radius, z);
0074   eta = -log(tan(theta / 2.));
0075   return make_pair(eta, phi);
0076 }
0077
0078 double
0079 PHG4Utils::get_eta(const double radius, const double z)
0080 {
0081   double eta;
0082   double theta;
0083   theta = atan2(radius, fabs(z));
0084   eta = -log(tan(theta / 2.));
0085   if (z < 0)
0086   {
0087     eta = -eta;
0088   }
0089   return eta;
0090 }
0091
0092 void PHG4Utils::SetColour(G4VisAttributes* att, const string& material)
0093 {
0094   if (!att)
0095   {
0096     cout << "G4VisAttributes pointer is NULL" << endl;
0097     return;
0098   }
0099   if (material == "AL_BABAR_MAG")
0100   {
0101     att->SetColour(G4Colour::Blue());
0102   }
0103   else if (material == "BlackHole" || material == "G4_AIR")
0104   {
0105     att->SetColour(G4Colour::Black());
0106   }
0107   else if (material == "C4F10")
0108   {
0109     att->SetColour(0., 0., 0.5, 0.25);
0110   }
0111   else if (material == "CF4")
0112   {
0113     att->SetColour(G4Colour::Magenta());
0114   }
0115   else if (material == "G4_Al")
0116   {
0117     att->SetColour(G4Colour::Gray());
0118   }
0119   else if (material == "G4_Au" || material == "G4_KAPTON" || material == "G4_Si")
0120   {
0121     att->SetColour(G4Colour::Yellow());
0122   }
0123   else if (material == "G4_CARBON_DIOXIDE" || material == "Quartz")
0124   {
0125     att->SetColour(G4Colour::Green());
0126   }
0127   else if (material == "G4_CELLULOSE_CELLOPHANE")
0128   {
0129     att->SetColour(0.25, 0.25, 0.);
0130   }
0131   else if (material == "G4_Cu")
0132   {
0133     att->SetColour(1., 0.51, 0.278);
0134   }
0135   else if (material == "G4_Fe")
0136   {
0137     att->SetColour(0.29, 0.44, 0.54);
0138   }
0139   else if (material == "G4_MYLAR")
0140   {
0141     att->SetColour(0.5, 0.5, 0.5, 0.25);
0142   }
0143   else if (material == "G4_METHANE")
0144   {
0145     att->SetColour(0., 1., 1., 0.25);
0146   }
0147   else if (material == "G4_TEFLON")
0148   {
0149     att->SetColour(G4Colour::White());
0150   }
0151   else if (material == "G4_W")
0152   {
0153     att->SetColour(0.36, 0.36, 0.36);
0154   }
0155   else if (material == "Scintillator" || material == "G4_POLYSTYRENE")
0156   {
0157     att->SetColour(0., 1., 1.);
0158   }
0159   else if (material == "W_Epoxy")
0160   {
0161     att->SetColour(0.5, 0.5, 0.5);
0162   }
0163   else if (material == "G10")
0164   {
0165     att->SetColour(1., 1., 0., 0.5);
0166   }
0167   else
0168   {
0169     // cout << "default color red for material " << material << endl;
0170     att->SetColour(G4Colour::Cyan());
0171   }
0172   return;
0173 }
0174 // returns success/failure and intersection point (x/y)
0175 std::pair<bool, std::pair<double, double>> PHG4Utils::lines_intersect(
0176     double ax,
0177     double ay,
0178     double bx,
0179     double by,
0180     double cx,
0181     double cy,
0182     double dx,
0183     double dy)
0184 {
0185   // Find if a line segment limited by points A and B
0186   // intersects line segment limited by points C and D.
0187   // First check if an infinite line defined by A and B intersects
0188   // segment (C,D). If h is from 0 to 1 line and line segment intersect
0189   // Then check in intersection point is between C and D
0190   double ex = bx - ax;  // E=B-A
0191   double ey = by - ay;
0192   double fx = dx - cx;  // F=D-C
0193   double fy = dy - cy;
0194   double px = -ey;  // P
0195   double py = ex;
0196
0197   double bottom = fx * px + fy * py;  // F*P
0198   double gx = ax - cx;                // A-C
0199   double gy = ay - cy;
0200   double top = gx * px + gy * py;  // G*P
0201
0202   double h = 99999.;
0203   if (bottom != 0.)
0204   {
0205     h = top / bottom;
0206   }
0207
0208   // intersection point R = C + F*h
0209   if (h > 0. && h < 1.)
0210   {
0211     double rx = cx + fx * h;
0212     double ry = cy + fy * h;
0213     // cout << "      line/segment intersection coordinates: " << *rx << " " << *ry << endl;
0214     if ((rx > ax && rx > bx) || (rx < ax && rx < bx) || (ry < ay && ry < by) || (ry > ay && ry > by))
0215     {
0216       // cout << "       NO segment/segment intersection!" << endl;
0217       return make_pair(false, make_pair(NAN, NAN));
0218     }
0219     else
0220     {
0221       // cout << "       segment/segment intersection!" << endl;
0222       return make_pair(true, make_pair(rx, ry));
0223     }
0224   }
0225
0226   return make_pair(false, make_pair(NAN, NAN));
0227 }
0228
0229 // returns success/failure and length of the line segment inside the rectangle (output)
0230 std::pair<bool, double> PHG4Utils::line_and_rectangle_intersect(
0231     double ax,
0232     double ay,
0233     double bx,
0234     double by,
0235     double cx,
0236     double cy,
0237     double dx,
0238     double dy)
0239 {
0240   // find if a line isegment limited by points (A,B)
0241   // intersects with a rectangle defined by two
0242   // corner points (C,D) two other points are E and F
0243   //   E--------D
0244   //   |        |
0245   //   |        |
0246   //   C--------F
0247
0248   if (cx > dx || cy > dy)
0249   {
0250     cout << PHWHERE << "ERROR: Bad rectangle definition!" << endl;
0251     return make_pair(false, NAN);
0252   }
0253
0254   double ex = cx;
0255   double ey = dy;
0256   double fx = dx;
0257   double fy = cy;
0258
0259   vector<double> vx;
0260   vector<double> vy;
0261   pair<bool, pair<double, double>> intersect1 = lines_intersect(ax, ay, bx, by, cx, cy, fx, fy);
0262   //  bool i1 = lines_intersect(ax, ay, bx, by, cx, cy, fx, fy, &rx, &ry);
0263   if (intersect1.first)
0264   {
0265     vx.push_back(intersect1.second.first);
0266     vy.push_back(intersect1.second.second);
0267   }
0268   pair<bool, pair<double, double>> intersect2 = lines_intersect(ax, ay, bx, by, fx, fy, dx, dy);
0269   //  bool i2 = lines_intersect(ax, ay, bx, by, fx, fy, dx, dy, &rx, &ry);
0270   if (intersect2.first)
0271   {
0272     vx.push_back(intersect2.second.first);
0273     vy.push_back(intersect2.second.second);
0274   }
0275   pair<bool, pair<double, double>> intersect3 = lines_intersect(ax, ay, bx, by, ex, ey, dx, dy);
0276   //  bool i3 = lines_intersect(ax, ay, bx, by, ex, ey, dx, dy, &rx, &ry);
0277   if (intersect3.first)
0278   {
0279     vx.push_back(intersect3.second.first);
0280     vy.push_back(intersect3.second.second);
0281   }
0282   pair<bool, pair<double, double>> intersect4 = lines_intersect(ax, ay, bx, by, cx, cy, ex, ey);
0283   //  bool i4 = lines_intersect(ax, ay, bx, by, cx, cy, ex, ey, &rx, &ry);
0284   if (intersect4.first)
0285   {
0286     vx.push_back(intersect4.second.first);
0287     vy.push_back(intersect4.second.second);
0288   }
0289
0290   // cout << "Rectangle intersections: " << i1 << " " << i2 << " " << i3 << " " << i4 << endl;
0291   // cout << "Number of intersections = " << vx.size() << endl;
0292
0293   double rr = 0.;
0294   if (vx.size() == 2)
0295   {
0296     rr = sqrt((vx[0] - vx[1]) * (vx[0] - vx[1]) + (vy[0] - vy[1]) * (vy[0] - vy[1]));
0297     //  cout << "Length of intersection = " << *rr << endl;
0298   }
0299   if (vx.size() == 1)
0300   {
0301     // find which point (A or B) is within the rectangle
0302     if (ax > cx && ay > cy && ax < dx && ay < dy)  // point A is inside the rectangle
0303     {
0304       // cout << "Point A is inside the rectangle." << endl;
0305       rr = sqrt((vx[0] - ax) * (vx[0] - ax) + (vy[0] - ay) * (vy[0] - ay));
0306     }
0307     if (bx > cx && by > cy && bx < dx && by < dy)  // point B is inside the rectangle
0308     {
0309       // cout << "Point B is inside the rectangle." << endl;
0310       rr = sqrt((vx[0] - bx) * (vx[0] - bx) + (vy[0] - by) * (vy[0] - by));
0311     }
0312   }
0313
0314   if (intersect1.first || intersect2.first || intersect3.first || intersect4.first)
0315   {
0316     return make_pair(true, rr);
0317   }
0318   return make_pair(false, NAN);
0319 }
0320
0321 // NOLINTNEXTLINE(misc-no-recursion)
0322 double PHG4Utils::sA(double r, double x, double y)
0323 {
0324   // Uses analytic formula for the integral of a circle between limits set by the corner of a rectangle
0325   // It is called repeatedly to find the overlap area between the circle and rectangle
0326   // I found this code implementing the integral on a web forum called "ars technica",
0327   // https://arstechnica.com/civis/viewtopic.php?t=306492
0328   // posted by "memp"
0329
0330   double a;
0331
0332   if (x < 0)
0333   {
0334     return -sA(r, -x, y);
0335   }
0336
0337   if (y < 0)
0338   {
0339     return -sA(r, x, -y);
0340   }
0341
0342   if (x > r)
0343   {
0344     x = r;
0345   }
0346
0347   if (y > r)
0348   {
0349     y = r;
0350   }
0351
0352   if (x * x + y * y > r * r)
0353   {
0354     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;
0355
0356     a *= 0.5;
0357   }
0358   else
0359   {
0360     a = x * y;
0361   }
0362
0363   return a;
0364 }
0365
0366 double PHG4Utils::circle_rectangle_intersection(double x1, double y1, double x2, double y2, double mx, double my, double r)
0367 {
0368   // Find the area of overlap of a circle and rectangle
0369   // Calls sA, which uses an analytic formula to determine the integral of the circle between limits set by the corners of the rectangle
0370
0371   // move the rectangle to the frame where the circle is at (0,0)
0372   x1 -= mx;
0373   x2 -= mx;
0374   y1 -= my;
0375   y2 -= my;
0376
0377   // {
0378   //   cout << " mx " << mx << " my " << my << " r " << r << " x1 " << x1 << " x2 " << x2 << " y1 " << y1 << " y2 " << y2 << endl;
0379   //   cout << " sA21 " << sA(r, x2, y1)
0380   //        << " sA11 " << sA(r, x1, y1)
0381   //        << " sA22 " << sA(r, x2, y2)
0382   //        << " sA12 " << sA(r, x1, y2)
0383   //        << endl;
0384   // }
0385
0386   return sA(r, x2, y1) - sA(r, x1, y1) - sA(r, x2, y2) + sA(r, x1, y2);
0387 }
0388
0389 string PHG4Utils::md5sum(const std::string& filename)
0390 {
0391   ifstream myfile;
0392   myfile.open(filename);
0393   if (!myfile.is_open())
0394   {
0395     cout << "Error opening " << filename << endl;
0396     gSystem->Exit(1);
0397     exit(1);
0398   }
0399   boost::uuids::detail::md5 hash;
0400   boost::uuids::detail::md5::digest_type digest;
0401   char c;
0402   while (myfile.get(c))
0403   {
0404     hash.process_bytes(&c, 1);
0405   }
0406   myfile.close();
0407   hash.get_digest(digest);
0408   const auto charDigest = reinterpret_cast<const char*>(&digest);
0409   std::string result;
0410   boost::algorithm::hex(charDigest, charDigest + sizeof(boost::uuids::detail::md5::digest_type), std::back_inserter(result));
0411   boost::algorithm::to_lower(result);
0412   return result;
0413 }