Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:19:00 

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