Back to home page

sPhenix code displayed by LXR
 
 

    

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

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