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