Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-17 09:21:39 

0001 #include "PHG4InnerHcalDetector.h"
0002 
0003 #include "PHG4HcalDefs.h"
0004 #include "PHG4InnerHcalDisplayAction.h"
0005 #include "PHG4InnerHcalSubsystem.h"
0006 
0007 #include <phparameter/PHParameters.h>
0008 
0009 #include <g4main/PHG4Detector.h>
0010 #include <g4main/PHG4DisplayAction.h>
0011 #include <g4main/PHG4Subsystem.h>
0012 #include <g4main/PHG4Utils.h>
0013 
0014 #include <phool/PHCompositeNode.h>
0015 #include <phool/PHIODataNode.h>
0016 #include <phool/PHNode.h>  // for PHNode
0017 #include <phool/PHNodeIterator.h>
0018 #include <phool/PHObject.h>  // for PHObject
0019 #include <phool/getClass.h>
0020 #include <phool/phool.h>
0021 #include <phool/recoConsts.h>
0022 
0023 #include <calobase/RawTowerDefs.h>           // for convert_name_...
0024 #include <calobase/RawTowerGeom.h>           // for RawTowerGeom
0025 #include <calobase/RawTowerGeomContainer.h>  // for RawTowerGeomC...
0026 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
0027 #include <calobase/RawTowerGeomv1.h>
0028 
0029 #include <TSystem.h>
0030 
0031 #include <Geant4/G4AssemblyVolume.hh>
0032 #include <Geant4/G4Box.hh>
0033 #include <Geant4/G4ExtrudedSolid.hh>
0034 #include <Geant4/G4IntersectionSolid.hh>
0035 #include <Geant4/G4LogicalVolume.hh>
0036 #include <Geant4/G4Material.hh>
0037 #include <Geant4/G4PVPlacement.hh>
0038 #include <Geant4/G4RotationMatrix.hh>
0039 #include <Geant4/G4String.hh>
0040 #include <Geant4/G4SystemOfUnits.hh>
0041 #include <Geant4/G4ThreeVector.hh>
0042 #include <Geant4/G4Transform3D.hh>
0043 #include <Geant4/G4Tubs.hh>
0044 #include <Geant4/G4TwoVector.hh>
0045 #include <Geant4/G4Types.hh>
0046 #include <Geant4/G4UserLimits.hh>
0047 #include <Geant4/G4VSolid.hh>
0048 
0049 #include <CGAL/Boolean_set_operations_2.h>
0050 #include <CGAL/Circular_kernel_intersections.h>
0051 #include <CGAL/Exact_circular_kernel_2.h>
0052 #include <CGAL/Object.h>
0053 #include <CGAL/point_generators_2.h>
0054 
0055 #include <boost/lexical_cast.hpp>
0056 #include <boost/tokenizer.hpp>
0057 
0058 #include <algorithm>
0059 #include <cmath>
0060 #include <cstdlib>
0061 #include <iostream>
0062 #include <iterator>
0063 #include <sstream>
0064 
0065 class PHCompositeNode;
0066 
0067 using Circle_2 = CGAL::Circle_2<PHG4InnerHcalDetector::Circular_k>;
0068 using Circular_arc_point_2 = CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>;
0069 using Line_2 = CGAL::Line_2<PHG4InnerHcalDetector::Circular_k>;
0070 using Segment_2 = CGAL::Segment_2<PHG4InnerHcalDetector::Circular_k>;
0071 #if CGAL_VERSION_NR > 1060000000
0072 typedef typename CGAL::CK2_Intersection_traits<PHG4InnerHcalDetector::Circular_k, Circle_2,Line_2>::type
0073        Intersection_result;
0074 #endif
0075 
0076 // there is still a minute problem for very low tilt angles where the scintillator
0077 // face touches the boundary instead of the corner, subtracting 1 permille from the total
0078 // scintilator length takes care of this
0079 namespace
0080 {
0081   double constexpr subtract_from_scinti_x = 0.1 * mm;
0082 }
0083 
0084 PHG4InnerHcalDetector::PHG4InnerHcalDetector(PHG4Subsystem *subsys, PHCompositeNode *Node, PHParameters *parameters, const std::string &dnam)
0085   : PHG4Detector(subsys, Node, dnam)
0086   , m_DisplayAction(dynamic_cast<PHG4InnerHcalDisplayAction *>(subsys->GetDisplayAction()))
0087   , m_Params(parameters)
0088   , m_InnerRadius(m_Params->get_double_param("inner_radius") * cm)
0089   , m_OuterRadius(m_Params->get_double_param("outer_radius") * cm)
0090   , m_SizeZ(m_Params->get_double_param("size_z") * cm)
0091   , m_ScintiTileZ(m_SizeZ)
0092   , m_ScintiTileThickness(m_Params->get_double_param("scinti_tile_thickness") * cm)
0093   , m_ScintiInnerGap(m_Params->get_double_param("scinti_inner_gap") * cm)
0094   , m_ScintiOuterGap(m_Params->get_double_param("scinti_outer_gap") * cm)
0095   , m_ScintiOuterRadius(m_Params->get_double_param("scinti_outer_radius") * cm)
0096   , m_TiltAngle(m_Params->get_double_param("tilt_angle") * deg)
0097   , m_EnvelopeInnerRadius(m_InnerRadius)
0098   , m_EnvelopeOuterRadius(m_OuterRadius)
0099   , m_EnvelopeZ(m_SizeZ)
0100   , m_NumScintiPlates(m_Params->get_int_param(PHG4HcalDefs::scipertwr) * m_Params->get_int_param("n_towers"))
0101   , m_NumScintiTilesPos(m_Params->get_int_param("n_scinti_tiles_pos"))
0102   , m_NumScintiTilesNeg(m_Params->get_int_param("n_scinti_tiles_neg"))
0103   , m_Active(m_Params->get_int_param("active"))
0104   , m_AbsorberActive(m_Params->get_int_param("absorberactive"))
0105   , m_ScintiLogicNamePrefix("HcalInnerScinti")
0106 {
0107   // n_scinti_tiles takes precedence
0108 
0109   int nTiles = 2 * m_Params->get_int_param(PHG4HcalDefs::n_scinti_tiles);
0110   if (nTiles <= 0)
0111   {
0112     nTiles = m_NumScintiTilesPos + m_NumScintiTilesNeg;
0113   }
0114   else
0115   {
0116     m_NumScintiTilesPos = nTiles / 2;
0117     m_Params->set_int_param("n_scinti_tiles_pos", nTiles / 2);
0118     m_NumScintiTilesNeg = nTiles / 2;
0119     m_Params->set_int_param("n_scinti_tiles_neg", nTiles / 2);
0120   }
0121 
0122   // allocate memory for scintillator plates
0123   m_ScintiTilesVec.assign(nTiles, static_cast<G4VSolid *>(nullptr));
0124 }
0125 
0126 PHG4InnerHcalDetector::~PHG4InnerHcalDetector()
0127 {
0128   delete m_ScintiMotherAssembly;
0129 }
0130 
0131 //_______________________________________________________________
0132 //_______________________________________________________________
0133 int PHG4InnerHcalDetector::IsInInnerHcal(G4VPhysicalVolume *volume) const
0134 {
0135   if (m_AbsorberActive)
0136   {
0137     if (m_SteelAbsorberPhysVolSet.contains(volume))
0138     {
0139       return -1;
0140     }
0141   }
0142   if (m_Active)
0143   {
0144     if (m_ScintiTilePhysVolMap.contains(volume))
0145     {
0146       return 1;
0147     }
0148   }
0149   return 0;
0150 }
0151 
0152 G4VSolid *
0153 PHG4InnerHcalDetector::ConstructScintillatorBox(G4LogicalVolume * /*hcalenvelope*/)
0154 {
0155   double mid_radius = m_InnerRadius + (m_OuterRadius - m_InnerRadius) / 2.;
0156   Point_2 p_in_1(mid_radius, 0);  // center of scintillator
0157 
0158   // length of upper edge (middle till outer circle intersect
0159   // x/y coordinate of end of center vertical
0160   double xcoord = m_ScintiTileThickness / 2. * sin(fabs(m_TiltAngle) / rad) + mid_radius;
0161   double ycoord = m_ScintiTileThickness / 2. * cos(fabs(m_TiltAngle) / rad) + 0;
0162   Point_2 p_upperedge(xcoord, ycoord);
0163   Line_2 s2(p_in_1, p_upperedge);  // center vertical
0164 
0165   Line_2 perp = s2.perpendicular(p_upperedge);
0166   Point_2 sc1(m_OuterRadius, 0);
0167   Point_2 sc2(0, m_OuterRadius);
0168   Point_2 sc3(-m_OuterRadius, 0);
0169   Circle_2 outer_circle(sc1, sc2, sc3);
0170 #if CGAL_VERSION_NR > 1060000000
0171     std::vector<Intersection_result> res;
0172 #else
0173   std::vector<CGAL::Object> res;
0174 #endif
0175   CGAL::intersection(outer_circle, perp, std::back_inserter(res));
0176   Point_2 upperright;
0177   for (const auto& obj : res)
0178   {
0179     if (const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned> *point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned>>(&obj))
0180     {
0181       if (CGAL::to_double(point->first.x()) > CGAL::to_double(p_upperedge.x()))
0182       {
0183         double deltax = CGAL::to_double(point->first.x()) - CGAL::to_double(p_upperedge.x());
0184         double deltay = CGAL::to_double(point->first.y()) - CGAL::to_double(p_upperedge.y());
0185         // the scintillator is twice as long
0186         m_ScintiTileXUpper = sqrt(deltax * deltax + deltay * deltay);  //
0187         Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
0188         upperright = pntmp;
0189       }
0190     }
0191     else
0192     {
0193       std::cout << "CGAL::Object type not pair..." << std::endl;
0194     }
0195   }
0196   // length of lower edge (middle till inner circle intersect
0197   xcoord = mid_radius - m_ScintiTileThickness / 2. * sin(fabs(m_TiltAngle) / rad);
0198   ycoord = 0 - m_ScintiTileThickness / 2. * cos(fabs(m_TiltAngle) / rad);
0199   Point_2 p_loweredge(xcoord, ycoord);
0200   Line_2 s3(p_in_1, p_loweredge);
0201   Line_2 l_lower = s3.perpendicular(p_loweredge);
0202   Point_2 ic1(m_InnerRadius, 0);
0203   Point_2 ic2(0, m_InnerRadius);
0204   Point_2 ic3(-m_InnerRadius, 0);
0205   Circle_2 inner_circle(ic1, ic2, ic3);
0206   res.clear();
0207   CGAL::intersection(inner_circle, l_lower, std::back_inserter(res));
0208   Point_2 lowerleft;
0209   // we have 2 intersections - we want the one furthest to the right (largest x). The correct one is
0210   // certainly > 0 but the second one depends on the tilt angle and might also be > 0
0211   double minx = 0;
0212   for (const auto& obj : res)
0213   {
0214     if (const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned> *point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned>>(&obj))
0215     {
0216       if (CGAL::to_double(point->first.x()) > minx)
0217       {
0218         minx = CGAL::to_double(point->first.x());
0219         double deltax = CGAL::to_double(point->first.x()) - CGAL::to_double(p_loweredge.x());
0220         double deltay = CGAL::to_double(point->first.y()) - CGAL::to_double(p_loweredge.y());
0221         m_ScintiTileXLower = sqrt(deltax * deltax + deltay * deltay);
0222         Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
0223         lowerleft = pntmp;
0224       }
0225     }
0226   }
0227   m_ScintiTileX = m_ScintiTileXUpper + m_ScintiTileXLower - ((m_OuterRadius - m_ScintiOuterRadius) / cos(m_TiltAngle / rad));
0228   m_ScintiTileX -= subtract_from_scinti_x;
0229   G4VSolid *scintibox = new G4Box("ScintiTile", m_ScintiTileX / 2., m_ScintiTileThickness / 2., m_ScintiTileZ / 2.);
0230   m_VolumeScintillator = scintibox->GetCubicVolume() * m_NumScintiPlates;
0231   return scintibox;
0232 }
0233 
0234 G4VSolid *
0235 PHG4InnerHcalDetector::ConstructSteelPlate(G4LogicalVolume * /*hcalenvelope*/)
0236 {
0237   // calculate steel plate on top of the scinti box. Lower edge is the upper edge of
0238   // the scintibox + 1/2 the airgap
0239   double mid_radius = m_InnerRadius + (m_OuterRadius - m_InnerRadius) / 2.;
0240   // first the lower edge, just like the scinti box, just add the air gap
0241   // and calculate intersection of edge with inner and outer radius.
0242   Point_2 p_in_1(mid_radius, 0);  // center of lower scintillator
0243   double angle_mid_scinti = M_PI / 2. + m_TiltAngle / rad;
0244   double xcoord = m_ScintiInnerGap / 2. * cos(angle_mid_scinti / rad) + mid_radius;
0245   double ycoord = m_ScintiInnerGap / 2. * sin(angle_mid_scinti / rad) + 0;
0246   Point_2 p_loweredge(xcoord, ycoord);
0247   Line_2 s2(p_in_1, p_loweredge);               // center vertical
0248   Line_2 perp = s2.perpendicular(p_loweredge);  // that is the lower edge of the steel plate
0249   Point_2 sc1(m_InnerRadius, 0);
0250   Point_2 sc2(0, m_InnerRadius);
0251   Point_2 sc3(-m_InnerRadius, 0);
0252   Circle_2 inner_circle(sc1, sc2, sc3);
0253 #if CGAL_VERSION_NR > 1060000000
0254     std::vector<Intersection_result> res;
0255 #else
0256   std::vector<CGAL::Object> res;
0257 #endif
0258   CGAL::intersection(inner_circle, perp, std::back_inserter(res));
0259   Point_2 lowerleft;
0260   for (const auto& obj : res)
0261   {
0262     if (const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned> *point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned>>(&obj))
0263     {
0264       if (CGAL::to_double(point->first.x()) > 0)
0265       {
0266         Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
0267         lowerleft = pntmp;
0268       }
0269     }
0270     else
0271     {
0272       std::cout << "CGAL::Object type not pair..." << std::endl;
0273     }
0274   }
0275 
0276   double xcoord2 = m_ScintiOuterGap / 2. * cos(angle_mid_scinti / rad) + mid_radius;
0277   double ycoord2 = m_ScintiOuterGap / 2. * sin(angle_mid_scinti / rad) + 0;
0278   Point_2 p_loweredge2(xcoord2, ycoord2);
0279   Line_2 s2_2(p_in_1, p_loweredge2);                // center vertical
0280   Line_2 perp2 = s2_2.perpendicular(p_loweredge2);  // that is the lower edge of the steel plate
0281 
0282   Point_2 so1(m_OuterRadius, 0);
0283   Point_2 so2(0, m_OuterRadius);
0284   Point_2 so3(-m_OuterRadius, 0);
0285   Circle_2 outer_circle(so1, so2, so3);
0286   res.clear();  // just clear the content from the last intersection search
0287   CGAL::intersection(outer_circle, perp2, std::back_inserter(res));
0288   Point_2 lowerright;
0289   for (const auto& obj : res)
0290   {
0291     if (const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned> *point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned>>(&obj))
0292     {
0293       if (CGAL::to_double(point->first.x()) > CGAL::to_double(p_loweredge.x()))
0294       {
0295         Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
0296         lowerright = pntmp;
0297       }
0298     }
0299     else
0300     {
0301       std::cout << "CGAL::Object type not pair..." << std::endl;
0302     }
0303   }
0304   // now we have the lower left and rigth corner, now find the upper edge
0305   // find the center of the upper scintilator
0306 
0307   double phi_midpoint = 2 * M_PI / m_NumScintiPlates;
0308   double xmidpoint = cos(phi_midpoint) * mid_radius;
0309   double ymidpoint = sin(phi_midpoint) * mid_radius;
0310   // angle of perp line at center of scintillator
0311   angle_mid_scinti = (M_PI / 2. - phi_midpoint) - (M_PI / 2. + m_TiltAngle / rad);
0312   double xcoordup = xmidpoint - m_ScintiInnerGap / 2. * sin(angle_mid_scinti / rad);
0313   double ycoordup = ymidpoint - m_ScintiInnerGap / 2. * cos(angle_mid_scinti / rad);
0314   Point_2 upperleft;
0315   Point_2 upperright;
0316   Point_2 mid_upperscint(xmidpoint, ymidpoint);
0317   Point_2 p_upperedge(xcoordup, ycoordup);
0318   Line_2 sup(mid_upperscint, p_upperedge);        // center vertical
0319   Line_2 perpA = sup.perpendicular(p_upperedge);  // that is the upper edge of the steel plate
0320   Point_2 sc1A(m_InnerRadius, 0);
0321   Point_2 sc2A(0, m_InnerRadius);
0322   Point_2 sc3A(-m_InnerRadius, 0);
0323   Circle_2 inner_circleA(sc1A, sc2A, sc3A);
0324 #if CGAL_VERSION_NR > 1060000000
0325     std::vector<Intersection_result> resA;
0326 #else
0327   std::vector<CGAL::Object> resA;
0328 #endif
0329   CGAL::intersection(inner_circleA, perpA, std::back_inserter(resA));
0330   double pxmax = 0.;
0331   for (const auto& obj : resA)
0332   {
0333     if (const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned> *point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned>>(&obj))
0334     {
0335       if (CGAL::to_double(point->first.x()) > pxmax)
0336       {
0337         pxmax = CGAL::to_double(point->first.x());
0338         Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
0339         upperleft = pntmp;
0340       }
0341     }
0342     else
0343     {
0344       std::cout << "CGAL::Object type not pair..." << std::endl;
0345     }
0346   }
0347 
0348   double xcoordup2 = xmidpoint - m_ScintiOuterGap / 2. * sin(angle_mid_scinti / rad);
0349   double ycoordup2 = ymidpoint - m_ScintiOuterGap / 2. * cos(angle_mid_scinti / rad);
0350   Point_2 p_upperedge2(xcoordup2, ycoordup2);
0351   Line_2 sup2(mid_upperscint, p_upperedge2);         // center vertical
0352   Line_2 perpA2 = sup2.perpendicular(p_upperedge2);  // that is the upper edge of the steel plate
0353 
0354   Point_2 so1A(m_OuterRadius, 0);
0355   Point_2 so2A(0, m_OuterRadius);
0356   Point_2 so3A(-m_OuterRadius, 0);
0357   Circle_2 outer_circleA(so1A, so2A, so3A);
0358   resA.clear();  // just clear the content from the last intersection search
0359   CGAL::intersection(outer_circleA, perpA2, std::back_inserter(resA));
0360   for (const auto& obj : resA)
0361   {
0362     if (const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned> *point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned>>(&obj))
0363     {
0364       if (CGAL::to_double(point->first.x()) > CGAL::to_double(p_loweredge.x()))
0365       {
0366         Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
0367         upperright = pntmp;
0368       }
0369     }
0370     else
0371     {
0372       std::cout << "CGAL::Object type not pair..." << std::endl;
0373     }
0374   }
0375   // the left corners are on a secant with the inner boundary, they need to be shifted
0376   // to be a tangent at the center
0377   ShiftSecantToTangent(lowerleft, upperleft, upperright, lowerright);
0378   G4TwoVector v1(CGAL::to_double(upperleft.x()), CGAL::to_double(upperleft.y()));
0379   G4TwoVector v2(CGAL::to_double(upperright.x()), CGAL::to_double(upperright.y()));
0380   G4TwoVector v3(CGAL::to_double(lowerright.x()), CGAL::to_double(lowerright.y()));
0381   G4TwoVector v4(CGAL::to_double(lowerleft.x()), CGAL::to_double(lowerleft.y()));
0382   std::vector<G4TwoVector> vertexes;
0383   vertexes.push_back(v1);
0384   vertexes.push_back(v2);
0385   vertexes.push_back(v3);
0386   vertexes.push_back(v4);
0387   G4TwoVector zero(0, 0);
0388   G4VSolid *steel_plate = new G4ExtrudedSolid("SteelPlate",
0389                                               vertexes,
0390                                               m_SizeZ / 2.0,
0391                                               zero, 1.0,
0392                                               zero, 1.0);
0393 
0394   //  DisplayVolume(steel_plate, hcalenvelope);
0395   m_VolumeSteel = steel_plate->GetCubicVolume() * m_NumScintiPlates;
0396   return steel_plate;
0397 }
0398 
0399 void PHG4InnerHcalDetector::ShiftSecantToTangent(Point_2 &lowleft, Point_2 &upleft, Point_2 &upright, Point_2 &lowright) const
0400 {
0401   Line_2 secant(lowleft, upleft);
0402   Segment_2 upedge(upleft, upright);
0403   Segment_2 lowedge(lowleft, lowright);
0404   double xmid = (CGAL::to_double(lowleft.x()) + CGAL::to_double(upleft.x())) / 2.;
0405   double ymid = (CGAL::to_double(lowleft.y()) + CGAL::to_double(upleft.y())) / 2.;
0406   Point_2 midpoint(xmid, ymid);
0407   Line_2 sekperp = secant.perpendicular(midpoint);
0408   Point_2 sc1(m_InnerRadius, 0);
0409   Point_2 sc2(0, m_InnerRadius);
0410   Point_2 sc3(-m_InnerRadius, 0);
0411   Circle_2 inner_circle(sc1, sc2, sc3);
0412 #if CGAL_VERSION_NR > 1060000000
0413     std::vector<Intersection_result> res;
0414 #else
0415   std::vector<CGAL::Object> res;
0416 #endif
0417   CGAL::intersection(inner_circle, sekperp, std::back_inserter(res));
0418   double pxmax = 0.;
0419   Point_2 tangtouch;
0420   for (const auto& obj : res)
0421   {
0422     if (const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned> *point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned>>(&obj))
0423     {
0424       if (CGAL::to_double(point->first.x()) > pxmax)
0425       {
0426         pxmax = CGAL::to_double(point->first.x());
0427         Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
0428         tangtouch = pntmp;
0429       }
0430     }
0431     else
0432     {
0433       std::cout << "CGAL::Object type not pair..." << std::endl;
0434     }
0435   }
0436   Line_2 leftside = sekperp.perpendicular(tangtouch);
0437   CGAL::Object result = CGAL::intersection(upedge, leftside);
0438   if (const Point_2 *ipoint = CGAL::object_cast<Point_2>(&result))
0439   {
0440     upleft = *ipoint;
0441   }
0442   result = CGAL::intersection(lowedge, leftside);
0443   if (const Point_2 *ipoint = CGAL::object_cast<Point_2>(&result))
0444   {
0445     lowleft = *ipoint;
0446   }
0447   return;
0448 }
0449 
0450 // Construct the envelope and the call the
0451 // actual inner hcal construction
0452 void PHG4InnerHcalDetector::ConstructMe(G4LogicalVolume *logicWorld)
0453 {
0454   recoConsts *rc = recoConsts::instance();
0455   G4Material *Air = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
0456   G4VSolid *hcal_envelope_cylinder = new G4Tubs("InnerHcal_envelope_solid", m_EnvelopeInnerRadius, m_EnvelopeOuterRadius, m_EnvelopeZ / 2., 0, 2 * M_PI);
0457   m_VolumeEnvelope = hcal_envelope_cylinder->GetCubicVolume();
0458   G4LogicalVolume *hcal_envelope_log = new G4LogicalVolume(hcal_envelope_cylinder, Air, G4String("Hcal_envelope"), nullptr, nullptr, nullptr);
0459 
0460   G4RotationMatrix hcal_rotm;
0461   hcal_rotm.rotateX(m_Params->get_double_param("rot_x") * deg);
0462   hcal_rotm.rotateY(m_Params->get_double_param("rot_y") * deg);
0463   hcal_rotm.rotateZ(m_Params->get_double_param("rot_z") * deg);
0464   G4VPhysicalVolume *mothervol = new G4PVPlacement(G4Transform3D(hcal_rotm, G4ThreeVector(m_Params->get_double_param("place_x") * cm, m_Params->get_double_param("place_y") * cm, m_Params->get_double_param("place_z") * cm)), hcal_envelope_log, "InnerHcalEnvelope", logicWorld, false, false, OverlapCheck());
0465   m_DisplayAction->SetMyTopVolume(mothervol);
0466   ConstructInnerHcal(hcal_envelope_log);
0467   std::vector<G4VPhysicalVolume *>::iterator it = m_ScintiMotherAssembly->GetVolumesIterator();
0468   for (unsigned int i = 0; i < m_ScintiMotherAssembly->TotalImprintedVolumes(); i++)
0469   {
0470     // G4AssemblyVolumes naming convention:
0471     //     av_WWW_impr_XXX_YYY_ZZZ
0472     // where:
0473 
0474     //     WWW - assembly volume instance number
0475     //     XXX - assembly volume imprint number
0476     //     YYY - the name of the placed logical volume
0477     //     ZZZ - the logical volume index inside the assembly volume
0478     // e.g. av_1_impr_82_HcalInnerScinti_11_pv_11
0479     // 82 the number of the scintillator mother volume
0480     // HcalInnerScinti_11: name of scintillator slat
0481     // 11: number of scintillator slat logical volume
0482     // use boost tokenizer to separate the _, then take value
0483     // after "impr" for mother volume and after "pv" for scintillator slat
0484     // use boost lexical cast for string -> int conversion
0485     // the CopyNo is the mother volume + scinti id
0486     // so we can use the CopyNo rather than decoding the string further
0487     // looking for "pv"
0488     boost::char_separator<char> sep("_");
0489     boost::tokenizer<boost::char_separator<char>> tok((*it)->GetName(), sep);
0490     boost::tokenizer<boost::char_separator<char>>::const_iterator tokeniter;
0491     for (tokeniter = tok.begin(); tokeniter != tok.end(); ++tokeniter)
0492     {
0493       if (*tokeniter == "impr")
0494       {
0495         ++tokeniter;
0496         if (tokeniter != tok.end())
0497         {
0498           int layer_id = boost::lexical_cast<int>(*tokeniter);
0499           // check detector description, for assemblyvolumes it is not possible
0500           // to give the first volume id=0, so they go from id=1 to id=n.
0501           // I am not going to start with fortran again - our indices start
0502           // at zero, id=0 to id=n-1. So subtract one here
0503           int tower_id = (*it)->GetCopyNo() - layer_id;
0504           layer_id--;
0505           std::pair<int, int> layer_twr = std::make_pair(layer_id, tower_id);
0506           m_ScintiTilePhysVolMap.insert(std::pair<G4VPhysicalVolume *, std::pair<int, int>>(*it, layer_twr));
0507           if (layer_id < 0 || layer_id >= m_NumScintiPlates)
0508           {
0509             std::cout << "invalid scintillator row " << layer_id
0510                       << ", valid range 0 < row < " << m_NumScintiPlates << std::endl;
0511             gSystem->Exit(1);
0512           }
0513         }
0514         else
0515         {
0516           std::cout << PHWHERE << " Error parsing " << (*it)->GetName()
0517                     << " for mother volume number " << std::endl;
0518           gSystem->Exit(1);
0519         }
0520         break;
0521       }
0522     }
0523     ++it;
0524   }
0525   if (!m_Params->get_int_param("saveg4hit"))
0526   {
0527     AddGeometryNode();
0528   }
0529   return;
0530 }
0531 
0532 int PHG4InnerHcalDetector::ConstructInnerHcal(G4LogicalVolume *hcalenvelope)
0533 {
0534   ConsistencyCheck();
0535   SetTiltViaNcross();  // if number of crossings is set, use it to determine tilt
0536   CheckTiltAngle();    // die if the tilt angle is out of range
0537   G4VSolid *steel_plate = ConstructSteelPlate(hcalenvelope);
0538   G4LogicalVolume *steel_logical = new G4LogicalVolume(steel_plate, GetDetectorMaterial(m_Params->get_string_param("material")), "HcalInnerSteelPlate", nullptr, nullptr, nullptr);
0539   m_DisplayAction->AddSteelVolume(steel_logical);
0540   m_ScintiMotherAssembly = ConstructHcalScintillatorAssembly(hcalenvelope);
0541   double phi = 0;
0542   double deltaphi = 2 * M_PI / m_NumScintiPlates;
0543   std::ostringstream name;
0544   double middlerad = m_OuterRadius - (m_OuterRadius - m_InnerRadius) / 2.;
0545   // for shorter scintillators we have to shift by 1/2 the "missing length"
0546   double scinti_tile_orig_length = m_ScintiTileXUpper + m_ScintiTileXLower - subtract_from_scinti_x;
0547   double shiftslat = fabs(m_ScintiTileXLower - m_ScintiTileXUpper) / 2. + (scinti_tile_orig_length - m_ScintiTileX) / 2.;
0548   // calculate phi offset (copied from code inside following loop):
0549   // first get the center point (phi=0) so it's middlerad/0
0550   // then shift the scintillator center as documented in loop
0551   // then
0552   // for positive tilt angles we need the lower left corner of the scintillator
0553   // for negative tilt angles we need the upper right corner of the scintillator
0554   // as it turns out the code uses the middle of the face of the scintillator
0555   // as reference, if this is a problem the code needs to be modified to
0556   // actually calculate the corner (but the math of the construction is that
0557   // the middle of the scintillator sits at zero)
0558   double xp = cos(phi) * middlerad;
0559   double yp = sin(phi) * middlerad;
0560   xp -= cos((-m_TiltAngle) / rad - phi) * shiftslat;
0561   yp += sin((-m_TiltAngle) / rad - phi) * shiftslat;
0562   if (m_TiltAngle > 0)
0563   {
0564     double xo = xp - (m_ScintiTileX / 2.) * cos(m_TiltAngle / rad);
0565     double yo = yp - (m_ScintiTileX / 2.) * sin(m_TiltAngle / rad);
0566     phi = -atan(yo / xo);
0567   }
0568   else if (m_TiltAngle < 0)
0569   {
0570     double xo = xp + (m_ScintiTileX / 2.) * cos(m_TiltAngle / rad);
0571     double yo = yp + (m_ScintiTileX / 2.) * sin(m_TiltAngle / rad);
0572     phi = -atan(yo / xo);
0573   }
0574   // else (for m_TiltAngle = 0) phi stays zero
0575   for (int i = 0; i < m_NumScintiPlates; i++)
0576   {
0577     G4RotationMatrix *Rot = new G4RotationMatrix();
0578     double ypos = sin(phi) * middlerad;
0579     double xpos = cos(phi) * middlerad;
0580     // the center of the scintillator is not the center of the inner hcal
0581     // but depends on the tilt angle. Therefore we need to shift
0582     // the center from the mid point
0583     ypos += sin((-m_TiltAngle) / rad - phi) * shiftslat;
0584     xpos -= cos((-m_TiltAngle) / rad - phi) * shiftslat;
0585     Rot->rotateZ(phi * rad + m_TiltAngle);
0586     G4ThreeVector g4vec(xpos, ypos, 0);
0587     // great the MakeImprint always adds 1 to the copy number and 0 has a
0588     // special meaning (which then also adds 1). Basically our volume names
0589     // will start at 1 instead of 0 and there is nothing short of patching this
0590     // method. I'll take care of this in the decoding of the volume name
0591     // AAAAAAARHGS
0592     m_ScintiMotherAssembly->MakeImprint(hcalenvelope, g4vec, Rot, i, OverlapCheck());
0593     delete Rot;
0594     Rot = new G4RotationMatrix();
0595     Rot->rotateZ(-phi * rad);
0596     name.str("");
0597     name << "InnerHcalSteel_" << i;
0598     m_SteelAbsorberPhysVolSet.insert(new G4PVPlacement(Rot, G4ThreeVector(0, 0, 0), steel_logical, name.str(), hcalenvelope, false, i, OverlapCheck()));
0599     phi += deltaphi;
0600   }
0601   return 0;
0602 }
0603 
0604 // split the big scintillator into tiles covering eta ranges
0605 // since they are tilted it is not a straightforward theta cut
0606 // it starts at a given eta at the inner radius but the outer radius needs adjusting
0607 void PHG4InnerHcalDetector::ConstructHcalSingleScintillators(G4LogicalVolume *hcalenvelope)
0608 {
0609   G4VSolid *bigtile = ConstructScintillatorBox(hcalenvelope);
0610 
0611   // eta->theta
0612   // scinti_eta_coverage takes precedence
0613   G4double eta_cov = m_Params->get_double_param("scinti_eta_coverage");
0614   if (eta_cov > 0)
0615   {
0616     m_Params->set_double_param("scinti_eta_coverage_pos", eta_cov);
0617     m_Params->set_double_param("scinti_eta_coverage_neg", eta_cov);
0618   }
0619   G4double delta_eta_pos = m_Params->get_double_param("scinti_eta_coverage_pos") / m_NumScintiTilesPos;
0620   G4double delta_eta_neg = m_Params->get_double_param("scinti_eta_coverage_neg") / m_NumScintiTilesNeg;
0621 
0622   G4double eta = 0;
0623   G4double theta;
0624   G4double x[4];
0625   G4double z[4];
0626   std::ostringstream name;
0627   double overhang = (m_ScintiTileX - (m_OuterRadius - m_InnerRadius)) / 2.;
0628   double offset = 1 * cm + overhang;  // add 1cm to make sure the G4ExtrudedSolid
0629   // is larger than the tile so we do not have
0630   // funny edge effects when overlapping vols
0631   double scinti_gap_neighbor = m_Params->get_double_param("scinti_gap_neighbor") * cm;
0632 
0633   // Positive first, then negative
0634 
0635   for (int i = 0; i < m_NumScintiTilesPos; i++)
0636   {
0637     theta = M_PI / 2 - PHG4Utils::get_theta(eta);  // theta = 90 for eta=0
0638     x[0] = m_InnerRadius - overhang;
0639     z[0] = tan(theta) * m_InnerRadius;
0640     x[1] = m_OuterRadius + overhang;  // since the tile is tilted, x is not at the outer radius but beyond
0641     z[1] = tan(theta) * m_OuterRadius;
0642     eta += delta_eta_pos;
0643     theta = M_PI / 2 - PHG4Utils::get_theta(eta);  // theta = 90 for eta=0
0644     x[2] = m_InnerRadius - overhang;
0645     z[2] = tan(theta) * m_InnerRadius;
0646     x[3] = m_OuterRadius + overhang;  // since the tile is tilted, x is not at the outer radius but beyond
0647     z[3] = tan(theta) * m_OuterRadius;
0648     // apply gap between scintillators
0649     z[0] += scinti_gap_neighbor / 2.;
0650     z[1] += scinti_gap_neighbor / 2.;
0651     z[2] -= scinti_gap_neighbor / 2.;
0652     z[3] -= scinti_gap_neighbor / 2.;
0653     Point_2 leftsidelow(z[0], x[0]);
0654     Point_2 leftsidehigh(z[1], x[1]);
0655     x[0] = m_InnerRadius - offset;
0656     z[0] = x_at_y(leftsidelow, leftsidehigh, x[0]);
0657     x[1] = m_OuterRadius + offset;
0658     z[1] = x_at_y(leftsidelow, leftsidehigh, x[1]);
0659     Point_2 rightsidelow(z[2], x[2]);
0660     Point_2 rightsidehigh(z[3], x[3]);
0661     x[2] = m_OuterRadius + offset;
0662     z[2] = x_at_y(rightsidelow, rightsidehigh, x[2]);
0663     x[3] = m_InnerRadius - offset;
0664     z[3] = x_at_y(rightsidelow, rightsidehigh, x[3]);
0665 
0666     std::vector<G4TwoVector> vertexes;
0667     for (int j = 0; j < 4; j++)
0668     {
0669       G4TwoVector v(x[j], z[j]);
0670       vertexes.push_back(v);
0671     }
0672     G4TwoVector zero(0, 0);
0673 
0674     G4VSolid *scinti = new G4ExtrudedSolid("ScintillatorTile",
0675                                            vertexes,
0676                                            m_ScintiTileThickness + 0.2 * mm,
0677                                            zero, 1.0,
0678                                            zero, 1.0);
0679     G4RotationMatrix *rotm = new G4RotationMatrix();
0680     rotm->rotateX(-90 * deg);
0681     name.str("");
0682     name << "scintillator_" << i << "_left";
0683     G4VSolid *scinti_tile = new G4IntersectionSolid(name.str(), bigtile, scinti, rotm, G4ThreeVector(-(m_InnerRadius + m_OuterRadius) / 2., 0, -m_Params->get_double_param("place_z") * cm));
0684     delete rotm;
0685     m_ScintiTilesVec[i + m_NumScintiTilesNeg] = scinti_tile;
0686   }
0687 
0688   eta = 0.0;  // reset
0689   for (int i = 0; i < m_NumScintiTilesNeg; i++)
0690   {
0691     theta = M_PI / 2 - PHG4Utils::get_theta(eta);  // theta = 90 for eta=0
0692     x[0] = m_InnerRadius - overhang;
0693     z[0] = tan(theta) * m_InnerRadius;
0694     x[1] = m_OuterRadius + overhang;  // since the tile is tilted, x is not at the outer radius but beyond
0695     z[1] = tan(theta) * m_OuterRadius;
0696     eta += delta_eta_neg;
0697     theta = M_PI / 2 - PHG4Utils::get_theta(eta);  // theta = 90 for eta=0
0698     x[2] = m_InnerRadius - overhang;
0699     z[2] = tan(theta) * m_InnerRadius;
0700     x[3] = m_OuterRadius + overhang;  // since the tile is tilted, x is not at the outer radius but beyond
0701     z[3] = tan(theta) * m_OuterRadius;
0702     // apply gap between scintillators
0703     z[0] += scinti_gap_neighbor / 2.;
0704     z[1] += scinti_gap_neighbor / 2.;
0705     z[2] -= scinti_gap_neighbor / 2.;
0706     z[3] -= scinti_gap_neighbor / 2.;
0707     Point_2 leftsidelow(z[0], x[0]);
0708     Point_2 leftsidehigh(z[1], x[1]);
0709     x[0] = m_InnerRadius - offset;
0710     z[0] = x_at_y(leftsidelow, leftsidehigh, x[0]);
0711     x[1] = m_OuterRadius + offset;
0712     z[1] = x_at_y(leftsidelow, leftsidehigh, x[1]);
0713     Point_2 rightsidelow(z[2], x[2]);
0714     Point_2 rightsidehigh(z[3], x[3]);
0715     x[2] = m_OuterRadius + offset;
0716     z[2] = x_at_y(rightsidelow, rightsidehigh, x[2]);
0717     x[3] = m_InnerRadius - offset;
0718     z[3] = x_at_y(rightsidelow, rightsidehigh, x[3]);
0719 
0720     std::vector<G4TwoVector> vertexes;
0721     for (int j = 0; j < 4; j++)
0722     {
0723       G4TwoVector v(x[j], z[j]);
0724       vertexes.push_back(v);
0725     }
0726     G4TwoVector zero(0, 0);
0727 
0728     G4VSolid *scinti = new G4ExtrudedSolid("ScintillatorTile",
0729                                            vertexes,
0730                                            m_ScintiTileThickness + 0.2 * mm,
0731                                            zero, 1.0,
0732                                            zero, 1.0);
0733 
0734     G4RotationMatrix *rotm = new G4RotationMatrix();
0735     rotm->rotateX(90 * deg);
0736     name.str("");
0737     name << "scintillator_" << i << "_right";
0738     G4VSolid *scinti_tile = new G4IntersectionSolid(name.str(), bigtile, scinti, rotm, G4ThreeVector(-(m_InnerRadius + m_OuterRadius) / 2., 0, -m_Params->get_double_param("place_z") * cm));
0739     m_ScintiTilesVec[m_NumScintiTilesNeg - i - 1] = scinti_tile;
0740     delete rotm;
0741   }
0742 
0743   // for (unsigned int i=0; i<m_ScintiTilesVec.size(); i++)
0744   //     {
0745   //       if (m_ScintiTilesVec[i])
0746   //     {
0747   //           DisplayVolume(m_ScintiTilesVec[i],hcalenvelope );
0748   //         }
0749   //     }
0750 
0751   return;
0752 }
0753 
0754 double
0755 PHG4InnerHcalDetector::x_at_y(Point_2 &p0, Point_2 &p1, double yin)
0756 {
0757   double xret = NAN;
0758   double x[2];
0759   double y[2];
0760   x[0] = CGAL::to_double(p0.x());
0761   y[0] = CGAL::to_double(p0.y());
0762   x[1] = CGAL::to_double(p1.x());
0763   y[1] = CGAL::to_double(p1.y());
0764   Line_2 l(p0, p1);
0765   double newx = fabs(x[0]) + fabs(x[1]);
0766   Point_2 p0new(-newx, yin);
0767   Point_2 p1new(newx, yin);
0768   Segment_2 seg(p0new, p1new);
0769   CGAL::Object result = CGAL::intersection(l, seg);
0770   if (const Point_2 *ipoint = CGAL::object_cast<Point_2>(&result))
0771   {
0772     xret = CGAL::to_double(ipoint->x());
0773   }
0774   else
0775   {
0776     std::cout << PHWHERE << " failed for y = " << y << std::endl;
0777     std::cout << "p0(x): " << CGAL::to_double(p0.x()) << ", p0(y): " << CGAL::to_double(p0.y()) << std::endl;
0778     std::cout << "p1(x): " << CGAL::to_double(p1.x()) << ", p1(y): " << CGAL::to_double(p1.y()) << std::endl;
0779     exit(1);
0780   }
0781   return xret;
0782 }
0783 
0784 G4AssemblyVolume *
0785 PHG4InnerHcalDetector::ConstructHcalScintillatorAssembly(G4LogicalVolume *hcalenvelope)
0786 {
0787   ConstructHcalSingleScintillators(hcalenvelope);
0788   G4AssemblyVolume *assmeblyvol = new G4AssemblyVolume();
0789   std::ostringstream name;
0790   G4ThreeVector g4vec;
0791 
0792   double steplimits = m_Params->get_double_param("steplimits") * cm;
0793   for (unsigned int i = 0; i < m_ScintiTilesVec.size(); i++)
0794   {
0795     name.str("");
0796     name << m_ScintiLogicNamePrefix << i;
0797     G4UserLimits *g4userlimits = nullptr;
0798     if (isfinite(steplimits))
0799     {
0800       g4userlimits = new G4UserLimits(steplimits);
0801     }
0802     G4LogicalVolume *scinti_tile_logic = new G4LogicalVolume(m_ScintiTilesVec[i], GetDetectorMaterial("G4_POLYSTYRENE"), name.str(), nullptr, nullptr, g4userlimits);
0803     m_DisplayAction->AddScintiVolume(scinti_tile_logic);
0804     assmeblyvol->AddPlacedVolume(scinti_tile_logic, g4vec, nullptr);
0805   }
0806   return assmeblyvol;
0807 }
0808 
0809 // check if tilt angle is reasonable - too large, no intersections with inner radius
0810 int PHG4InnerHcalDetector::CheckTiltAngle() const
0811 {
0812   if (fabs(m_TiltAngle / rad) >= M_PI)
0813   {
0814     std::cout << PHWHERE << "invalid tilt angle, abs(tilt) >= 90 deg: " << (m_TiltAngle / deg)
0815               << std::endl;
0816     exit(1);
0817   }
0818 
0819   double mid_radius = m_InnerRadius + (m_OuterRadius - m_InnerRadius) / 2.;
0820   Point_2 pmid(mid_radius, 0);  // center of scintillator
0821   double xcoord = 0;
0822   double ycoord = mid_radius * tan(m_TiltAngle / rad);
0823   Point_2 pxnull(xcoord, ycoord);
0824   Line_2 s2(pmid, pxnull);
0825   Point_2 sc1(m_InnerRadius, 0);
0826   Point_2 sc2(0, m_InnerRadius);
0827   Point_2 sc3(-m_InnerRadius, 0);
0828   Circle_2 inner_circle(sc1, sc2, sc3);
0829 #if CGAL_VERSION_NR > 1060000000
0830     std::vector<Intersection_result> res;
0831 #else
0832   std::vector<CGAL::Object> res;
0833 #endif
0834   CGAL::intersection(inner_circle, s2, std::back_inserter(res));
0835   if (res.empty())
0836   {
0837     std::cout << PHWHERE << " Tilt angle " << (m_TiltAngle / deg)
0838               << " too large, no intersection with inner radius" << std::endl;
0839     exit(1);
0840   }
0841   return 0;
0842 }
0843 
0844 int PHG4InnerHcalDetector::ConsistencyCheck() const
0845 {
0846   // just make sure the parameters make a bit of sense
0847   if (m_InnerRadius >= m_OuterRadius)
0848   {
0849     std::cout << PHWHERE << ": Inner Radius " << m_InnerRadius / cm
0850               << " cm larger than Outer Radius " << m_OuterRadius / cm
0851               << " cm" << std::endl;
0852     gSystem->Exit(1);
0853   }
0854   if (m_ScintiTileThickness > m_ScintiInnerGap)
0855   {
0856     std::cout << PHWHERE << "Scintillator thickness " << m_ScintiTileThickness / cm
0857               << " cm larger than scintillator inner gap " << m_ScintiInnerGap / cm
0858               << " cm" << std::endl;
0859     gSystem->Exit(1);
0860   }
0861   if (m_ScintiOuterRadius <= m_InnerRadius)
0862   {
0863     std::cout << PHWHERE << "Scintillator outer radius " << m_ScintiOuterRadius / cm
0864               << " cm smaller than inner radius " << m_InnerRadius / cm
0865               << " cm" << std::endl;
0866     gSystem->Exit(1);
0867   }
0868   return 0;
0869 }
0870 
0871 void PHG4InnerHcalDetector::SetTiltViaNcross()
0872 {
0873   // if tilt angle is set it takes precedence
0874   // this happens during readback where the tilt angle was set
0875   // in this method
0876   int ncross = m_Params->get_int_param("ncross");
0877   if (!ncross || isfinite(m_TiltAngle))
0878   {
0879     // flag ncrossing as not used
0880     m_Params->set_int_param("ncross", 0);
0881     return;
0882   }
0883   if ((isfinite(m_TiltAngle)) && (Verbosity() > 0))
0884   {
0885     std::cout << "both number of crossings and tilt angle are set" << std::endl;
0886     std::cout << "using number of crossings to determine tilt angle" << std::endl;
0887   }
0888   double mid_radius = m_InnerRadius + (m_OuterRadius - m_InnerRadius) / 2.;
0889   double deltaphi = (2 * M_PI / m_NumScintiPlates) * ncross;
0890   Point_2 pnull(0, 0);
0891   Point_2 plow(m_InnerRadius, 0);
0892   Point_2 phightmp(1, tan(deltaphi));
0893   Point_2 pin1(m_InnerRadius, 0);
0894   Point_2 pin2(0, m_InnerRadius);
0895   Point_2 pin3(-m_InnerRadius, 0);
0896   Circle_2 inner_circle(pin1, pin2, pin3);
0897   Point_2 pmid1(mid_radius, 0);
0898   Point_2 pmid2(0, mid_radius);
0899   Point_2 pmid3(-mid_radius, 0);
0900   Circle_2 mid_circle(pmid1, pmid2, pmid3);
0901   Point_2 pout1(m_OuterRadius, 0);
0902   Point_2 pout2(0, m_OuterRadius);
0903   Point_2 pout3(-m_OuterRadius, 0);
0904   Circle_2 outer_circle(pout1, pout2, pout3);
0905   Line_2 l_up(pnull, phightmp);
0906 #if CGAL_VERSION_NR > 1060000000
0907     std::vector<Intersection_result> res;
0908 #else
0909   std::vector<CGAL::Object> res;
0910 #endif
0911   CGAL::intersection(outer_circle, l_up, std::back_inserter(res));
0912   Point_2 upperright;
0913   for (const auto& obj : res)
0914   {
0915     if (const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned> *point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned>>(&obj))
0916     {
0917       if (CGAL::to_double(point->first.x()) > 0)
0918       {
0919         Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
0920         upperright = pntmp;
0921       }
0922     }
0923     else
0924     {
0925       std::cout << "CGAL::Object type not pair..." << std::endl;
0926       exit(1);
0927     }
0928   }
0929   Line_2 l_right(plow, upperright);
0930   res.clear();
0931   Point_2 midpoint;
0932   CGAL::intersection(mid_circle, l_right, std::back_inserter(res));
0933   for (const auto& obj : res)
0934   {
0935     if (const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned> *point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned>>(&obj))
0936     {
0937       if (CGAL::to_double(point->first.x()) > 0)
0938       {
0939         Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
0940         midpoint = pntmp;
0941       }
0942     }
0943     else
0944     {
0945       std::cout << "CGAL::Object type not pair..." << std::endl;
0946       exit(1);
0947     }
0948   }
0949   // length left side
0950   double ll = sqrt((CGAL::to_double(midpoint.x()) - m_InnerRadius) * (CGAL::to_double(midpoint.x()) - m_InnerRadius) + CGAL::to_double(midpoint.y()) * CGAL::to_double(midpoint.y()));
0951   double upside = sqrt(CGAL::to_double(midpoint.x()) * CGAL::to_double(midpoint.x()) + CGAL::to_double(midpoint.y()) * CGAL::to_double(midpoint.y()));
0952   //  c^2 = a^2+b^2 - 2ab*cos(gamma)
0953   // gamma = acos((a^2+b^2=c^2)/2ab
0954   double tiltangle = acos((ll * ll + upside * upside - m_InnerRadius * m_InnerRadius) / (2 * ll * upside));
0955   tiltangle = tiltangle * rad;
0956   m_TiltAngle = copysign(tiltangle, ncross);
0957   m_Params->set_double_param("tilt_angle", m_TiltAngle / deg);
0958   return;
0959 }
0960 
0961 void PHG4InnerHcalDetector::Print(const std::string &what) const
0962 {
0963   std::cout << "Inner Hcal Detector:" << std::endl;
0964   if (what == "ALL" || what == "VOLUME")
0965   {
0966     std::cout << "Volume Envelope: " << m_VolumeEnvelope / cm3 << " cm^3" << std::endl;
0967     std::cout << "Volume Steel: " << m_VolumeSteel / cm3 << " cm^3" << std::endl;
0968     std::cout << "Volume Scintillator: " << m_VolumeScintillator / cm3 << " cm^3" << std::endl;
0969     std::cout << "Volume Air: " << (m_VolumeEnvelope - m_VolumeSteel - m_VolumeScintillator) / cm3 << " cm^3" << std::endl;
0970   }
0971   return;
0972 }
0973 
0974 std::pair<int, int> PHG4InnerHcalDetector::GetLayerTowerId(G4VPhysicalVolume *volume) const
0975 {
0976   auto it = m_ScintiTilePhysVolMap.find(volume);
0977   if (it != m_ScintiTilePhysVolMap.end())
0978   {
0979     return it->second;
0980   }
0981   std::cout << "could not locate volume " << volume->GetName()
0982             << " in Inner Hcal scintillator map" << std::endl;
0983   gSystem->Exit(1);
0984   // that's dumb but code checkers do not know that gSystem->Exit()
0985   // terminates, so using the standard exit() makes them happy
0986   exit(1);
0987 }
0988 
0989 // This is dulplicated code, we can get rid of it when we have the code to make towergeom for real data reco.
0990 void PHG4InnerHcalDetector::AddGeometryNode()
0991 {
0992   PHNodeIterator iter(topNode());
0993   PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0994   if (!runNode)
0995   {
0996     std::cout << PHWHERE << "Run Node missing, exiting." << std::endl;
0997     gSystem->Exit(1);
0998     exit(1);
0999   }
1000   PHNodeIterator runIter(runNode);
1001   PHCompositeNode *RunDetNode = dynamic_cast<PHCompositeNode *>(runIter.findFirst("PHCompositeNode", m_SuperDetector));
1002   if (!RunDetNode)
1003   {
1004     RunDetNode = new PHCompositeNode(m_SuperDetector);
1005     runNode->addNode(RunDetNode);
1006   }
1007   m_TowerGeomNodeName = "TOWERGEOM_" + m_SuperDetector;
1008   m_RawTowerGeom = findNode::getClass<RawTowerGeomContainer>(topNode(), m_TowerGeomNodeName);
1009   if (!m_RawTowerGeom)
1010   {
1011     m_RawTowerGeom = new RawTowerGeomContainer_Cylinderv1(RawTowerDefs::convert_name_to_caloid(m_SuperDetector));
1012     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(m_RawTowerGeom, m_TowerGeomNodeName, "PHObject");
1013     RunDetNode->addNode(newNode);
1014   }
1015   double innerrad = m_Params->get_double_param(PHG4HcalDefs::innerrad);
1016   double thickness = m_Params->get_double_param(PHG4HcalDefs::outerrad) - innerrad;
1017   m_RawTowerGeom->set_radius(innerrad);
1018   m_RawTowerGeom->set_thickness(thickness);
1019   m_RawTowerGeom->set_phibins(m_Params->get_int_param(PHG4HcalDefs::n_towers));
1020   m_RawTowerGeom->set_etabins(m_Params->get_int_param("etabins"));
1021   double geom_ref_radius = innerrad + thickness / 2.;
1022   double phistart = m_Params->get_double_param("phistart");
1023   if (!std::isfinite(phistart))
1024   {
1025     std::cout << PHWHERE << " phistart is not finite: " << phistart
1026               << ", exiting now (this will crash anyway)" << std::endl;
1027     gSystem->Exit(1);
1028   }
1029   for (int i = 0; i < m_Params->get_int_param(PHG4HcalDefs::n_towers); i++)
1030   {
1031     double phiend = phistart + 2. * M_PI / m_Params->get_int_param(PHG4HcalDefs::n_towers);
1032     std::pair<double, double> range = std::make_pair(phiend, phistart);
1033     phistart = phiend;
1034     m_RawTowerGeom->set_phibounds(i, range);
1035   }
1036   double etalowbound = -m_Params->get_double_param("scinti_eta_coverage_neg");
1037   for (int i = 0; i < m_Params->get_int_param("etabins"); i++)
1038   {
1039     // double etahibound = etalowbound + 2.2 / get_int_param("etabins");
1040     double etahibound = etalowbound +
1041                         (m_Params->get_double_param("scinti_eta_coverage_neg") + m_Params->get_double_param("scinti_eta_coverage_pos")) / m_Params->get_int_param("etabins");
1042     std::pair<double, double> range = std::make_pair(etalowbound, etahibound);
1043     m_RawTowerGeom->set_etabounds(i, range);
1044     etalowbound = etahibound;
1045   }
1046   for (int iphi = 0; iphi < m_RawTowerGeom->get_phibins(); iphi++)
1047   {
1048     for (int ieta = 0; ieta < m_RawTowerGeom->get_etabins(); ieta++)
1049     {
1050       const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::convert_name_to_caloid(m_SuperDetector), ieta, iphi);
1051 
1052       const double x(geom_ref_radius * cos(m_RawTowerGeom->get_phicenter(iphi)));
1053       const double y(geom_ref_radius * sin(m_RawTowerGeom->get_phicenter(iphi)));
1054       const double z(geom_ref_radius / tan(PHG4Utils::get_theta(m_RawTowerGeom->get_etacenter(ieta))));
1055 
1056       RawTowerGeom *tg = m_RawTowerGeom->get_tower_geometry(key);
1057       if (tg)
1058       {
1059         if (Verbosity() > 0)
1060         {
1061           std::cout << "IHCalDetector::InitRun - Tower geometry " << key << " already exists" << std::endl;
1062         }
1063 
1064         if (fabs(tg->get_center_x() - x) > 1e-4)
1065         {
1066           std::cout << "IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key << " with existing x = " << tg->get_center_x() << " and expected x = " << x
1067                     << std::endl;
1068 
1069           return;
1070         }
1071         if (fabs(tg->get_center_y() - y) > 1e-4)
1072         {
1073           std::cout << "IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key << " with existing y = " << tg->get_center_y() << " and expected y = " << y
1074                     << std::endl;
1075           return;
1076         }
1077         if (fabs(tg->get_center_z() - z) > 1e-4)
1078         {
1079           std::cout << "IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key << " with existing z= " << tg->get_center_z() << " and expected z = " << z
1080                     << std::endl;
1081           return;
1082         }
1083       }
1084       else
1085       {
1086         if (Verbosity() > 0)
1087         {
1088           std::cout << "IHCalDetector::InitRun - building tower geometry " << key << "" << std::endl;
1089         }
1090 
1091         tg = new RawTowerGeomv1(key);
1092 
1093         tg->set_center_x(x);
1094         tg->set_center_y(y);
1095         tg->set_center_z(z);
1096         m_RawTowerGeom->add_tower_geometry(tg);
1097       }
1098     }
1099   }
1100   if (Verbosity() > 0)
1101   {
1102     m_RawTowerGeom->identify();
1103   }
1104 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!