Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:17:50 

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