Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 // $$Id: PHG4FullProjTiltedSpacalDetector.cc,v 1.3 2015/02/10 15:39:26 pinkenbu Exp $$
0002 
0003 /*!
0004  * \file ${file_name}
0005  * \brief
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $$Revision: 1.3 $$
0008  * \date $$Date: 2015/02/10 15:39:26 $$
0009  */
0010 #include "PHG4FullProjTiltedSpacalDetector.h"
0011 
0012 #include "PHG4SpacalDisplayAction.h"
0013 
0014 #include <g4gdml/PHG4GDMLConfig.hh>
0015 
0016 #include <phool/recoConsts.h>
0017 
0018 #include <phparameter/PHParameters.h>
0019 
0020 #include <Geant4/G4Box.hh>
0021 #include <Geant4/G4DisplacedSolid.hh>
0022 #include <Geant4/G4Exception.hh>          // for G4Exception
0023 #include <Geant4/G4ExceptionSeverity.hh>  // for FatalException
0024 #include <Geant4/G4LogicalVolume.hh>
0025 #include <Geant4/G4PVPlacement.hh>
0026 #include <Geant4/G4PhysicalConstants.hh>
0027 #include <Geant4/G4String.hh>  // for G4String
0028 #include <Geant4/G4SystemOfUnits.hh>
0029 #include <Geant4/G4ThreeVector.hh>  // for G4ThreeVector
0030 #include <Geant4/G4Transform3D.hh>  // for G4Transform3D, G4TranslateY3D
0031 #include <Geant4/G4Trap.hh>
0032 #include <Geant4/G4Types.hh>  // for G4double
0033 #include <Geant4/G4Vector3D.hh>
0034 
0035 #include <TSystem.h>
0036 
0037 #include <algorithm>
0038 #include <cassert>
0039 #include <cmath>
0040 #include <iostream>  // for operator<<, basic_ostream
0041 #include <limits>    // for numeric_limits
0042 #include <map>       // for map<>::value_type, map
0043 #include <memory>    // for allocator_traits<>::value_...
0044 #include <numeric>   // std::accumulate
0045 #include <sstream>
0046 #include <string>  // std::string, std::to_string
0047 #include <vector>  // for vector
0048 
0049 class G4Material;
0050 class G4VSolid;
0051 class PHCompositeNode;
0052 
0053 //_______________________________________________________________
0054 // note this inactive thickness is ~1.5% of a radiation length
0055 PHG4FullProjTiltedSpacalDetector::PHG4FullProjTiltedSpacalDetector(PHG4Subsystem* subsys, PHCompositeNode* Node,
0056                                                                    const std::string& dnam, PHParameters* parameters, const int lyr)
0057   : PHG4SpacalDetector(subsys, Node, dnam, parameters, lyr, false)
0058   , m_Params(parameters)
0059 {
0060   assert(_geom == nullptr);
0061 
0062   _geom = new SpacalGeom_t();
0063   if (_geom == nullptr)
0064   {
0065     std::cout
0066         << "PHG4FullProjTiltedSpacalDetector::Constructor - Fatal Error - invalid geometry object!"
0067         << std::endl;
0068     gSystem->Exit(1);
0069   }
0070   assert(parameters);
0071 
0072   assert(parameters);
0073   get_geom_v3()->ImportParameters(*parameters);
0074   //  std::cout <<"PHG4FullProjTiltedSpacalDetector::Constructor -  get_geom_v3()->Print();"<<std::endl;
0075   //  get_geom_v3()->Print();
0076 }
0077 
0078 //_______________________________________________________________
0079 void PHG4FullProjTiltedSpacalDetector::ConstructMe(G4LogicalVolume* logicWorld)
0080 {
0081   if (get_geom_v3()->get_construction_verbose() >= 1)
0082   {
0083     std::cout << "PHG4FullProjTiltedSpacalDetector::Construct::" << GetName()
0084               << " - start with PHG4SpacalDetector::Construct()." << std::endl;
0085   }
0086 
0087   PHG4SpacalDetector::ConstructMe(logicWorld);
0088 
0089   if (get_geom_v3()->get_construction_verbose() >= 1)
0090   {
0091     std::cout << "PHG4FullProjTiltedSpacalDetector::Construct::" << GetName()
0092               << " - Completed." << std::endl;
0093   }
0094 
0095   if (m_Params->get_int_param("saveg4hit"))
0096   {
0097     return;
0098   }
0099   try
0100   {
0101     AddCellGeometryNode();
0102     AddTowerGeometryNode();
0103   }
0104   catch (std::exception& e)
0105   {
0106     std::cout << e.what() << std::endl;
0107     // exit(1);
0108   }
0109 }
0110 
0111 std::pair<G4LogicalVolume*, G4Transform3D>
0112 PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg()
0113 {
0114   if (!(get_geom_v3()->get_azimuthal_n_sec() > 4))
0115   {
0116     std::cout << "azimuthal n sec <= 4: " << get_geom_v3()->get_azimuthal_n_sec() << std::endl;
0117     gSystem->Exit(1);
0118   }
0119 
0120   // basic tilt geometry
0121   const G4double half_chord_backend =
0122       get_geom_v3()->get_max_radius() * cm * tan(pi / get_geom_v3()->get_azimuthal_n_sec())  //
0123       + fabs(get_geom_v3()->get_thickness() * cm * 0.5 * tan(get_geom_v3()->get_azimuthal_tilt()));
0124   const G4double reduced_outer_radius = sqrt(pow(get_geom_v3()->get_max_radius() * cm, 2) - half_chord_backend * half_chord_backend);
0125   const G4double enclosure_depth = reduced_outer_radius - get_geom_v3()->get_radius() * cm;
0126   const G4double enclosure_center = 0.5 * (reduced_outer_radius + get_geom_v3()->get_radius() * cm);
0127   const G4double enclosure_half_height_half_width = enclosure_center * tan(pi / get_geom_v3()->get_azimuthal_n_sec());
0128 
0129   const G4double width_adj1 = tan(get_geom_v3()->get_azimuthal_tilt() - pi / get_geom_v3()->get_azimuthal_n_sec()) * enclosure_depth * 0.5;
0130   const G4double width_adj2 = tan(get_geom_v3()->get_azimuthal_tilt() + pi / get_geom_v3()->get_azimuthal_n_sec()) * enclosure_depth * 0.5;
0131 
0132   const G4double center_adj = (width_adj1 + width_adj2) * 0.5;
0133   const G4double center_tilt_angle = atan2(center_adj, enclosure_depth * 0.5);
0134   const G4double inner_half_width = enclosure_half_height_half_width + 0.5 * (width_adj1 - width_adj2);
0135   const G4double outter_half_width = enclosure_half_height_half_width + 0.5 * (-width_adj1 + width_adj2);
0136 
0137   // enclosure walls
0138   const G4double edge1_tilt_angle = atan2(width_adj1, enclosure_depth * 0.5);
0139   const G4double edge2_tilt_angle = atan2(width_adj2, enclosure_depth * 0.5);
0140   const G4double edge1_half_depth = sqrt(width_adj1 * width_adj1 + enclosure_depth * enclosure_depth * .25);
0141   const G4double edge2_half_depth = sqrt(width_adj2 * width_adj2 + enclosure_depth * enclosure_depth * .25);
0142 
0143   // projective center
0144   const G4double half_projection_ratio = 0.5 * (-width_adj1 + width_adj2) / enclosure_half_height_half_width;
0145   const G4double projection_center_y = enclosure_center - ((enclosure_depth * 0.5) / half_projection_ratio);
0146   const G4double projection_center_x = center_adj / half_projection_ratio;
0147 
0148   // blocks azimuthal segmentation
0149   const int phi_bin_in_sec = get_geom_v3()->get_max_phi_bin_in_sec();
0150   assert(phi_bin_in_sec >= 1);
0151   const G4double block_azimuth_angle = (edge2_tilt_angle - edge1_tilt_angle) / phi_bin_in_sec;
0152   assert(block_azimuth_angle > 0);
0153   if (!(fabs(block_azimuth_angle - M_PI * 2 / get_geom_v3()->get_azimuthal_n_sec() / phi_bin_in_sec) < M_PI * std::numeric_limits<G4double>::epsilon()))
0154   {
0155     std::cout << "angle/nsec out of range: " << M_PI * std::numeric_limits<G4double>::epsilon() << std::endl;
0156     gSystem->Exit(1);
0157   }
0158   const G4double block_edge1_half_width = enclosure_half_height_half_width - (get_geom_v3()->get_sidewall_thickness() * cm + get_geom_v3()->get_sidewall_outer_torr() * cm + 2.0 * get_geom_v3()->get_assembly_spacing() * cm) / cos(edge1_tilt_angle);
0159   const G4double block_edge2_half_width = enclosure_half_height_half_width - (get_geom_v3()->get_sidewall_thickness() * cm + get_geom_v3()->get_sidewall_outer_torr() * cm + 2.0 * get_geom_v3()->get_assembly_spacing() * cm) / cos(edge2_tilt_angle);
0160   G4double block_width_ratio = 0;
0161   for (int sa = 0; sa < phi_bin_in_sec; ++sa)
0162   {
0163     block_width_ratio += 1 / cos(block_azimuth_angle * (0.5 + sa) + edge1_tilt_angle);
0164   }
0165   const G4double block_half_height_width = (block_edge1_half_width + block_edge2_half_width) / block_width_ratio;
0166   assert(block_half_height_width > 0);
0167 
0168   // write out the azimuthal block geometry
0169   // block azimuth geometry records
0170   struct block_azimuth_geom
0171   {
0172     G4double angle;
0173     G4double projection_center_y;
0174     G4double projection_center_x;
0175     G4double projection_length;
0176   };
0177   std::vector<block_azimuth_geom> block_azimuth_geoms(phi_bin_in_sec,
0178                                                       block_azimuth_geom{
0179                                                           std::numeric_limits<double>::quiet_NaN(),
0180                                                           std::numeric_limits<double>::quiet_NaN(),
0181                                                           std::numeric_limits<double>::quiet_NaN(),
0182                                                           std::numeric_limits<double>::quiet_NaN()});  // [phi-bin in sector] -> azimuth geometry
0183   G4double block_x_edge1 = block_edge1_half_width;
0184   for (int sa = 0; sa < phi_bin_in_sec; ++sa)
0185   {
0186     block_azimuth_geom& geom = block_azimuth_geoms[sa];
0187 
0188     geom.angle = block_azimuth_angle * (0.5 + sa) + edge1_tilt_angle;
0189     const G4double block_x_size = block_half_height_width / cos(geom.angle);
0190     assert(block_x_size > 0);
0191     const G4double x_center = block_x_edge1 - 0.5 * block_x_size;
0192 
0193     // projection center per block
0194     geom.projection_length = block_half_height_width / 2. / tan(block_azimuth_angle / 2.);
0195     assert(geom.projection_length > 0);
0196     geom.projection_center_y = enclosure_center - geom.projection_length * cos(geom.angle);
0197     geom.projection_center_x = x_center + geom.projection_length * sin(geom.angle);
0198 
0199     // next step
0200     block_x_edge1 -= block_x_size;
0201   }
0202 
0203   // write out the azimuthal block divider's geometry
0204   struct block_divider_azimuth_geom
0205   {
0206     G4double angle;  //! rotation angle
0207     G4double projection_center_y;
0208     G4double projection_center_x;
0209     G4double thickness;            // thickness in the approximate azimuth direction
0210     G4double radial_displacement;  //! displacement along the width direction, which is the radial direction if tilt = 0
0211     G4double width;                //! wdith along the approximate radial direction
0212   };
0213   assert(phi_bin_in_sec >= 1);
0214   std::vector<block_divider_azimuth_geom> divider_azimuth_geoms(phi_bin_in_sec - 1,
0215                                                                 block_divider_azimuth_geom{
0216                                                                     std::numeric_limits<double>::quiet_NaN(),
0217                                                                     std::numeric_limits<double>::quiet_NaN(),
0218                                                                     std::numeric_limits<double>::quiet_NaN(),
0219                                                                     std::numeric_limits<double>::quiet_NaN(),
0220                                                                     std::numeric_limits<double>::quiet_NaN(),
0221                                                                     std::numeric_limits<double>::quiet_NaN()});
0222 
0223   if (get_geom_v3()->get_sidewall_thickness() > 0)
0224   {
0225     for (int sa = 0; sa < phi_bin_in_sec - 1; ++sa)
0226     {
0227       block_divider_azimuth_geom& geom = divider_azimuth_geoms[sa];
0228 
0229       geom.angle = 0.5 * (block_azimuth_geoms[sa].angle + block_azimuth_geoms[sa + 1].angle);
0230       geom.projection_center_y = 0.5 * (block_azimuth_geoms[sa].projection_center_y + block_azimuth_geoms[sa + 1].projection_center_y);
0231       geom.projection_center_x = 0.5 * (block_azimuth_geoms[sa].projection_center_x + block_azimuth_geoms[sa + 1].projection_center_x);
0232       geom.radial_displacement = 0.5 * (block_azimuth_geoms[sa].projection_length + block_azimuth_geoms[sa + 1].projection_length);
0233 
0234       geom.thickness = 2.0 * get_geom_v3()->get_assembly_spacing() * cm * cos(block_azimuth_angle / 2.) - 2 * um;
0235       geom.width = get_geom_v3()->get_divider_width() * cm;
0236     }
0237   }
0238 
0239   if (fabs(block_x_edge1 - (-block_edge2_half_width)) > get_geom_v3()->get_assembly_spacing() * cm)
0240   {
0241     std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg - ERROR - " << std::endl
0242               << "\t block_x_edge1 = " << block_x_edge1 << std::endl
0243               << "\t block_edge2_half_width = " << block_edge2_half_width << std::endl
0244               << "\t fabs(block_x_edge1 - (-block_edge2_half_width)) = " << fabs(block_x_edge1 - (-block_edge2_half_width)) << std::endl
0245               << "\t get_geom_v3()->get_assembly_spacing() * cm = " << get_geom_v3()->get_assembly_spacing() * cm << std::endl;
0246   }
0247   if (!(fabs(block_x_edge1 - (-block_edge2_half_width)) < get_geom_v3()->get_assembly_spacing() * cm))  // closure check
0248   {
0249     std::cout << "closure check failed: " << fabs(block_x_edge1 - (-block_edge2_half_width)) << std::endl;
0250     gSystem->Exit(1);
0251   }
0252 
0253   if (Verbosity())
0254   {
0255     std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg - " << std::endl
0256               << "\t edge1_tilt_angle = " << edge1_tilt_angle << std::endl
0257               << "\t edge2_tilt_angle = " << edge2_tilt_angle << std::endl
0258               << "\t projection_center_y = " << projection_center_y << std::endl
0259               << "\t projection_center_x = " << projection_center_x << std::endl
0260               << "\t block_azimuth_angle = " << block_azimuth_angle << std::endl
0261               << "\t block_edge1_half_width = " << block_edge1_half_width << std::endl
0262               << "\t block_edge2_half_width = " << block_edge2_half_width << std::endl
0263               << "\t block_width_ratio = " << block_width_ratio << std::endl
0264               << "\t block_half_height_width = " << block_half_height_width << std::endl;
0265 
0266     for (int sa = 0; sa < phi_bin_in_sec; ++sa)
0267     {
0268       std::cout << "\t block[" << sa << "].angle = " << block_azimuth_geoms[sa].angle << std::endl;
0269       std::cout << "\t block[" << sa << "].projection_center_y = " << block_azimuth_geoms[sa].projection_center_y << std::endl;
0270       std::cout << "\t block[" << sa << "].projection_center_x = " << block_azimuth_geoms[sa].projection_center_x << std::endl;
0271     }
0272     for (int sa = 0; sa < phi_bin_in_sec - 1; ++sa)
0273     {
0274       std::cout << "\t divider[" << sa << "].angle = " << divider_azimuth_geoms[sa].angle << std::endl;
0275       std::cout << "\t divider[" << sa << "].projection_center_x = " << divider_azimuth_geoms[sa].projection_center_x << std::endl;
0276       std::cout << "\t divider[" << sa << "].projection_center_y = " << divider_azimuth_geoms[sa].projection_center_y << std::endl;
0277       std::cout << "\t divider[" << sa << "].radial_displacement = " << divider_azimuth_geoms[sa].radial_displacement << std::endl;
0278       std::cout << "\t divider[" << sa << "].thickness = " << divider_azimuth_geoms[sa].thickness << std::endl;
0279       std::cout << "\t divider[" << sa << "].width = " << divider_azimuth_geoms[sa].width << std::endl;
0280     }
0281   }
0282 
0283   assert(enclosure_depth > 10 * cm);
0284 
0285   G4VSolid* sec_solid = new G4Trap(
0286       /*const G4String& pName*/ G4String(GetName() + std::string("_sec_trap")),
0287       enclosure_depth * 0.5,                                                                              // G4double pDz,
0288       center_tilt_angle, halfpi,                                                                          // G4double pTheta, G4double pPhi,
0289       inner_half_width, get_geom_v3()->get_length() * cm / 2.0, get_geom_v3()->get_length() * cm / 2.0,   // G4double pDy1, G4double pDx1, G4double pDx2,
0290       0,                                                                                                  // G4double pAlp1,
0291       outter_half_width, get_geom_v3()->get_length() * cm / 2.0, get_geom_v3()->get_length() * cm / 2.0,  // G4double pDy2, G4double pDx3, G4double pDx4,
0292       0                                                                                                   // G4double pAlp2 //
0293   );
0294   G4Transform3D sec_solid_transform = G4TranslateY3D(enclosure_center) * G4RotateY3D(halfpi) * G4RotateX3D(-halfpi);
0295   G4VSolid* sec_solid_place = new G4DisplacedSolid(G4String(GetName() + std::string("_sec")), sec_solid, sec_solid_transform);
0296 
0297   recoConsts* rc = recoConsts::instance();
0298   G4Material* cylinder_mat = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
0299   assert(cylinder_mat);
0300 
0301   G4LogicalVolume* sec_logic = new G4LogicalVolume(sec_solid_place, cylinder_mat,
0302                                                    G4String(G4String(GetName() + std::string("_sec"))), nullptr, nullptr, nullptr);
0303 
0304   GetDisplayAction()->AddVolume(sec_logic, "Sector");
0305 
0306   // construct walls
0307 
0308   G4Material* wall_mat = GetDetectorMaterial(get_geom_v3()->get_sidewall_mat());
0309   assert(wall_mat);
0310 
0311   if (get_geom_v3()->get_sidewall_thickness() > 0)
0312   {
0313     // end walls
0314     if (get_geom_v3()->get_construction_verbose() >= 1)
0315     {
0316       std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg::" << GetName()
0317                 << " - construct end walls." << std::endl;
0318     }
0319     //        G4Tubs* wall_solid = new G4Tubs(G4String(GetName() + std::string("_EndWall")),
0320     //            get_geom_v3()->get_radius() * cm + get_geom_v3()->get_sidewall_outer_torr() * cm,
0321     //            get_geom_v3()->get_max_radius() * cm - get_geom_v3()->get_sidewall_outer_torr() * cm,
0322     //            get_geom_v3()->get_sidewall_thickness() * cm / 2.0,
0323     //            halfpi - pi / get_geom_v3()->get_azimuthal_n_sec(),
0324     //            twopi / get_geom_v3()->get_azimuthal_n_sec());
0325     const G4double side_wall_half_thickness = get_geom_v3()->get_sidewall_thickness() * cm / 2.0;
0326     G4VSolid* wall_solid = new G4Trap(G4String(GetName() + std::string("_EndWall_trap")),
0327                                       enclosure_depth * 0.5,                                                  // G4double pDz,
0328                                       center_tilt_angle, halfpi,                                              // G4double pTheta, G4double pPhi,
0329                                       inner_half_width, side_wall_half_thickness, side_wall_half_thickness,   // G4double pDy1, G4double pDx1, G4double pDx2,
0330                                       0,                                                                      // G4double pAlp1,
0331                                       outter_half_width, side_wall_half_thickness, side_wall_half_thickness,  // G4double pDy2, G4double pDx3, G4double pDx4,
0332                                       0                                                                       // G4double pAlp2 //
0333     );
0334     G4VSolid* wall_solid_place = new G4DisplacedSolid(G4String(GetName() + std::string("_EndWall")), wall_solid, sec_solid_transform);
0335 
0336     G4LogicalVolume* wall_logic = new G4LogicalVolume(wall_solid_place, wall_mat,
0337                                                       G4String(G4String(GetName() + std::string("_EndWall"))), nullptr, nullptr,
0338                                                       nullptr);
0339     GetDisplayAction()->AddVolume(wall_logic, "Wall");
0340 
0341     using z_locations_t = std::map<int, double>;
0342     z_locations_t z_locations;
0343     z_locations[000] = get_geom_v3()->get_sidewall_thickness() * cm / 2.0 + get_geom_v3()->get_assembly_spacing() * cm;
0344     z_locations[001] = get_geom_v3()->get_length() * cm / 2.0 - (get_geom_v3()->get_sidewall_thickness() * cm / 2.0 + get_geom_v3()->get_assembly_spacing() * cm);
0345     z_locations[100] = -(get_geom_v3()->get_sidewall_thickness() * cm / 2.0 + get_geom_v3()->get_assembly_spacing() * cm);
0346     z_locations[101] = -(get_geom_v3()->get_length() * cm / 2.0 - (get_geom_v3()->get_sidewall_thickness() * cm / 2.0 + get_geom_v3()->get_assembly_spacing() * cm));
0347 
0348     for (z_locations_t::value_type& val : z_locations)
0349     {
0350       if (get_geom_v3()->get_construction_verbose() >= 2)
0351       {
0352         std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg::"
0353                   << GetName() << " - constructed End Wall ID " << val.first
0354                   << " @ Z = " << val.second << std::endl;
0355       }
0356       G4Transform3D wall_trans = G4TranslateZ3D(val.second);
0357 
0358       G4PVPlacement* wall_phys = new G4PVPlacement(wall_trans, wall_logic,
0359                                                    G4String(GetName()) + G4String("_EndWall_") + std::to_string(val.first), sec_logic,
0360                                                    false, val.first, OverlapCheck());
0361 
0362       calo_vol[wall_phys] = val.first;
0363       assert(gdml_config);
0364       gdml_config->exclude_physical_vol(wall_phys);
0365     }
0366   }
0367   //
0368   if (get_geom_v3()->get_sidewall_thickness() > 0)
0369   {
0370     // side walls
0371     if (get_geom_v3()->get_construction_verbose() >= 1)
0372     {
0373       std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg::" << GetName()
0374                 << " - construct side walls." << std::endl;
0375     }
0376 
0377     using sign_t = std::map<int, std::pair<int, int>>;
0378     sign_t signs;
0379     signs[100] = std::make_pair(+1, +1);
0380     signs[101] = std::make_pair(+1, -1);
0381     signs[200] = std::make_pair(-1, +1);
0382     signs[201] = std::make_pair(-1, -1);
0383 
0384     for (sign_t::value_type& val : signs)
0385     {
0386       const int sign_z = val.second.first;
0387       const int sign_azimuth = val.second.second;
0388 
0389       if (get_geom_v3()->get_construction_verbose() >= 2)
0390       {
0391         std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg::"
0392                   << GetName() << " - constructed Side Wall ID " << val.first
0393                   << " with"
0394                   << " Shift X = "
0395                   << sign_azimuth * (get_geom_v3()->get_sidewall_thickness() * cm / 2.0 + get_geom_v3()->get_sidewall_outer_torr() * cm)
0396                   << " Rotation Z = "
0397                   << sign_azimuth * pi / get_geom_v3()->get_azimuthal_n_sec()
0398                   << " Shift Z = " << sign_z * (get_geom_v3()->get_length() * cm / 4)
0399                   << std::endl;
0400       }
0401       const G4double azimuth_roate = sign_azimuth > 0 ? edge1_tilt_angle : edge2_tilt_angle;
0402       const G4double edge_half_depth = -get_geom_v3()->get_sidewall_thickness() * cm - get_geom_v3()->get_sidewall_outer_torr() * cm + (sign_azimuth > 0 ? edge1_half_depth : edge2_half_depth);
0403 
0404       G4Box* wall_solid = new G4Box(G4String(GetName() + G4String("_SideWall_") + std::to_string(val.first)),
0405                                     get_geom_v3()->get_sidewall_thickness() * cm / 2.0,
0406                                     edge_half_depth,
0407                                     (get_geom_v3()->get_length() / 2. - 2 * (get_geom_v3()->get_sidewall_thickness() + 2. * get_geom_v3()->get_assembly_spacing())) * cm * .5);
0408 
0409       G4LogicalVolume* wall_logic = new G4LogicalVolume(wall_solid, wall_mat,
0410                                                         G4String(G4String(GetName() + G4String("_SideWall_") + std::to_string(val.first))), nullptr, nullptr,
0411                                                         nullptr);
0412       GetDisplayAction()->AddVolume(wall_logic, "Wall");
0413 
0414       const G4Transform3D wall_trans =
0415           G4TranslateZ3D(sign_z * (get_geom_v3()->get_length() * cm / 4)) *
0416           G4TranslateY3D(enclosure_center) *
0417           G4TranslateX3D(sign_azimuth * enclosure_half_height_half_width) *
0418           G4RotateZ3D(azimuth_roate) *
0419           G4TranslateX3D(-sign_azimuth * (get_geom_v3()->get_sidewall_thickness() * cm / 2.0 + get_geom_v3()->get_sidewall_outer_torr() * cm));
0420 
0421       G4PVPlacement* wall_phys = new G4PVPlacement(wall_trans, wall_logic,
0422                                                    G4String(GetName()) + G4String("_SideWall_") + std::to_string(val.first), sec_logic,
0423                                                    false, val.first, OverlapCheck());
0424 
0425       calo_vol[wall_phys] = val.first;
0426       assert(gdml_config);
0427       gdml_config->exclude_physical_vol(wall_phys);
0428     }
0429   }
0430 
0431   // construct dividers
0432   if (get_geom_v3()->get_divider_width() > 0)
0433   {
0434     if (get_geom_v3()->get_construction_verbose() >= 1)
0435     {
0436       std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg::" << GetName()
0437                 << " - construct dividers" << std::endl;
0438     }
0439 
0440     G4Material* divider_mat = GetDetectorMaterial(get_geom_v3()->get_divider_mat());
0441     assert(divider_mat);
0442 
0443     int ID = 300;
0444     for (const auto& geom : divider_azimuth_geoms)
0445     {
0446       G4Box* divider_solid = new G4Box(G4String(GetName() + G4String("_Divider_") + std::to_string(ID)),
0447                                        geom.thickness / 2.0,
0448                                        geom.width / 2.,
0449                                        (get_geom_v3()->get_length() / 2. - 2 * (get_geom_v3()->get_sidewall_thickness() + 2. * get_geom_v3()->get_assembly_spacing())) * cm * .5);
0450 
0451       G4LogicalVolume* wall_logic = new G4LogicalVolume(divider_solid, divider_mat,
0452                                                         G4String(G4String(GetName() + G4String("_Divider_") + std::to_string(ID))), nullptr, nullptr,
0453                                                         nullptr);
0454       GetDisplayAction()->AddVolume(wall_logic, "Divider");
0455 
0456       for (int sign_z = -1; sign_z <= 1; sign_z += 2)
0457       {
0458         G4Transform3D wall_trans =
0459             G4TranslateX3D(geom.projection_center_x) *
0460             G4TranslateY3D(geom.projection_center_y) *
0461             G4RotateZ3D(geom.angle) *
0462             G4TranslateY3D(geom.radial_displacement) *
0463             G4TranslateZ3D(sign_z * (get_geom_v3()->get_length() * cm / 4));
0464 
0465         G4PVPlacement* wall_phys = new G4PVPlacement(wall_trans, wall_logic,
0466                                                      G4String(GetName()) + G4String("_Divider_") + std::to_string(ID), sec_logic,
0467                                                      false, ID, OverlapCheck());
0468 
0469         calo_vol[wall_phys] = ID;
0470         assert(gdml_config);
0471         gdml_config->exclude_physical_vol(wall_phys);
0472 
0473         if (Verbosity())
0474         {
0475           std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg - placing divider " << wall_phys->GetName() << " copy ID " << ID << std::endl;
0476         }
0477 
0478         ++ID;
0479       }
0480     }  //    for (const auto & geom : divider_azimuth_geoms)
0481   }
0482 
0483   //  // construct towers
0484   //
0485   for (const SpacalGeom_t::tower_map_t::value_type& val : get_geom_v3()->get_sector_tower_map())
0486   {
0487     SpacalGeom_t::geom_tower g_tower = val.second;
0488 
0489     const int tower_id = g_tower.id;
0490     const int tower_phi_id_in_sec = tower_id % 10;
0491     assert(tower_phi_id_in_sec >= 0);
0492     assert(tower_phi_id_in_sec < phi_bin_in_sec);
0493 
0494     const auto& block_azimuth_geom = block_azimuth_geoms.at(tower_phi_id_in_sec);
0495 
0496     G4LogicalVolume* LV_tower = Construct_Tower(g_tower);
0497 
0498     G4Transform3D block_trans =
0499         G4TranslateX3D(block_azimuth_geom.projection_center_x) *
0500         G4TranslateY3D(block_azimuth_geom.projection_center_y) *
0501         G4RotateZ3D(block_azimuth_geom.angle) *
0502         G4TranslateX3D(g_tower.centralX * cm) *
0503         G4TranslateY3D(g_tower.centralY * cm) *
0504         G4TranslateZ3D(g_tower.centralZ * cm) *
0505         G4RotateX3D(g_tower.pRotationAngleX * rad);
0506 
0507     const bool overlapcheck_block = OverlapCheck() and (get_geom_v3()->get_construction_verbose() >= 2);
0508 
0509     G4PVPlacement* block_phys = new G4PVPlacement(block_trans, LV_tower,
0510                                                   G4String(GetName()) + G4String("_Tower_") + std::to_string(g_tower.id), sec_logic, false,
0511                                                   g_tower.id, overlapcheck_block);
0512     block_vol[block_phys] = g_tower.id;
0513     assert(gdml_config);
0514     gdml_config->exclude_physical_vol(block_phys);
0515 
0516     if (g_tower.LightguideHeight > 0)
0517     {
0518       // also build a light guide
0519 
0520       for (int ix = 0; ix < g_tower.NSubtowerX; ix++)
0521       //  int ix = 0;
0522       {
0523         for (int iy = 0; iy < g_tower.NSubtowerY; iy++)
0524         //        int iy = 0;
0525         {
0526           G4LogicalVolume* LV_lg = Construct_LightGuide(g_tower, ix,
0527                                                         iy);
0528 
0529           G4PVPlacement* lg_phys = new G4PVPlacement(block_trans, LV_lg, LV_lg->GetName(),
0530                                                      sec_logic, false, g_tower.id, overlapcheck_block);
0531 
0532           block_vol[lg_phys] = g_tower.id * 100 + ix * 10 + iy;
0533 
0534           assert(gdml_config);
0535           gdml_config->exclude_physical_vol(lg_phys);
0536         }
0537       }
0538     }
0539   }
0540 
0541   std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg::" << GetName()
0542             << " - constructed " << get_geom_v3()->get_sector_tower_map().size()
0543             << " unique towers" << std::endl;
0544 
0545   return std::make_pair(sec_logic, G4Transform3D::Identity);
0546 }
0547 
0548 //! Fully projective spacal with 2D tapered modules. To speed up construction, same-length fiber is used cross one tower
0549 int PHG4FullProjTiltedSpacalDetector::Construct_Fibers_SameLengthFiberPerTower(
0550     const PHG4FullProjTiltedSpacalDetector::SpacalGeom_t::geom_tower& g_tower,
0551     G4LogicalVolume* LV_tower)
0552 {
0553   // construct fibers
0554 
0555   // first check out the fibers geometry
0556 
0557   using fiber_par_map = std::map<int, std::pair<G4Vector3D, G4Vector3D>>;
0558   fiber_par_map fiber_par;
0559   G4double min_fiber_length = g_tower.pDz * cm * 4;
0560 
0561   G4Vector3D v_zshift = G4Vector3D(tan(g_tower.pTheta) * cos(g_tower.pPhi),
0562                                    tan(g_tower.pTheta) * sin(g_tower.pPhi), 1) *
0563                         g_tower.pDz;
0564   //  int fiber_ID = 0;
0565   for (int ix = 0; ix < g_tower.NFiberX; ix++)
0566   //  int ix = 0;
0567   {
0568     const double weighted_ix = static_cast<double>(ix) / (g_tower.NFiberX - 1.);
0569 
0570     const double weighted_pDx1 = (g_tower.pDx1 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0571     const double weighted_pDx2 = (g_tower.pDx2 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0572 
0573     const double weighted_pDx3 = (g_tower.pDx3 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0574     const double weighted_pDx4 = (g_tower.pDx4 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0575 
0576     for (int iy = 0; iy < g_tower.NFiberY; iy++)
0577     //        int iy = 0;
0578     {
0579       if ((ix + iy) % 2 == 1)
0580       {
0581         continue;  // make a triangle pattern
0582       }
0583 
0584       const double weighted_iy = static_cast<double>(iy) / (g_tower.NFiberY - 1.);
0585 
0586       const double weighted_pDy1 = (g_tower.pDy1 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
0587       const double weighted_pDy2 = (g_tower.pDy2 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
0588 
0589       const double weighted_pDx12 = weighted_pDx1 * (1 - weighted_iy) + weighted_pDx2 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp1);
0590       const double weighted_pDx34 = weighted_pDx3 * (1 - weighted_iy) + weighted_pDx4 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp2);
0591 
0592       G4Vector3D v1 = G4Vector3D(weighted_pDx12, weighted_pDy1, 0) - v_zshift;
0593       G4Vector3D v2 = G4Vector3D(weighted_pDx34, weighted_pDy2, 0) + v_zshift;
0594 
0595       G4Vector3D vector_fiber = (v2 - v1);
0596       vector_fiber *= (vector_fiber.mag() - get_geom_v3()->get_fiber_outer_r()) / vector_fiber.mag();  // shrink by fiber boundary protection
0597       G4Vector3D center_fiber = (v2 + v1) / 2;
0598 
0599       // convert to Geant4 units
0600       vector_fiber *= cm;
0601       center_fiber *= cm;
0602 
0603       const int fiber_ID = g_tower.compose_fiber_id(ix, iy);
0604       fiber_par[fiber_ID] = std::make_pair(vector_fiber,
0605                                            center_fiber);
0606 
0607       const G4double fiber_length = vector_fiber.mag();
0608 
0609       min_fiber_length = std::min(fiber_length, min_fiber_length);
0610 
0611       //          ++fiber_ID;
0612     }
0613   }
0614 
0615   int fiber_count = 0;
0616 
0617   const G4double fiber_length = min_fiber_length;
0618   std::vector<G4double> fiber_cut;
0619 
0620   std::stringstream ss;
0621   ss << std::string("_Tower") << g_tower.id;
0622   G4LogicalVolume* fiber_logic = Construct_Fiber(fiber_length, ss.str());
0623 
0624   for (const fiber_par_map::value_type& val : fiber_par)
0625   {
0626     const int fiber_ID = val.first;
0627     G4Vector3D vector_fiber = val.second.first;
0628     G4Vector3D center_fiber = val.second.second;
0629     const G4double optimal_fiber_length = vector_fiber.mag();
0630 
0631     const G4Vector3D v1 = center_fiber - 0.5 * vector_fiber;
0632 
0633     // keep a statistics
0634     assert(optimal_fiber_length - fiber_length >= 0);
0635     fiber_cut.push_back(optimal_fiber_length - fiber_length);
0636 
0637     center_fiber += (fiber_length / optimal_fiber_length - 1) * 0.5 * vector_fiber;
0638     vector_fiber *= fiber_length / optimal_fiber_length;
0639 
0640     //      const G4Vector3D v1_new = center_fiber - 0.5 *vector_fiber;
0641 
0642     if (get_geom_v3()->get_construction_verbose() >= 3)
0643     {
0644       std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_Fibers_SameLengthFiberPerTower::" << GetName()
0645                 << " - constructed fiber " << fiber_ID << ss.str()  //
0646                 << ", Length = " << optimal_fiber_length << "-"
0647                 << (optimal_fiber_length - fiber_length) << "mm, "  //
0648                 << "x = " << center_fiber.x() << "mm, "             //
0649                 << "y = " << center_fiber.y() << "mm, "             //
0650                 << "z = " << center_fiber.z() << "mm, "             //
0651                 << "vx = " << vector_fiber.x() << "mm, "            //
0652                 << "vy = " << vector_fiber.y() << "mm, "            //
0653                 << "vz = " << vector_fiber.z() << "mm, "            //
0654                 << std::endl;
0655     }
0656 
0657     const G4double rotation_angle = G4Vector3D(0, 0, 1).angle(vector_fiber);
0658     const G4Vector3D rotation_axis =
0659         rotation_angle == 0 ? G4Vector3D(1, 0, 0) : G4Vector3D(0, 0, 1).cross(vector_fiber);
0660 
0661     G4Transform3D fiber_place(
0662         G4Translate3D(center_fiber.x(), center_fiber.y(), center_fiber.z()) * G4Rotate3D(rotation_angle, rotation_axis));
0663 
0664     std::stringstream name;
0665     name << GetName() + std::string("_Tower") << g_tower.id << "_fiber"
0666          << ss.str();
0667 
0668     const bool overlapcheck_fiber = OverlapCheck() and (get_geom_v3()->get_construction_verbose() >= 3);
0669     G4PVPlacement* fiber_physi = new G4PVPlacement(fiber_place, fiber_logic,
0670                                                    G4String(name.str()), LV_tower, false, fiber_ID,
0671                                                    overlapcheck_fiber);
0672     fiber_vol[fiber_physi] = fiber_ID;
0673     assert(gdml_config);
0674     gdml_config->exclude_physical_vol(fiber_physi);
0675 
0676     fiber_count++;
0677   }
0678 
0679   if (get_geom_v3()->get_construction_verbose() >= 2)
0680   {
0681     std::cout
0682         << "PHG4FullProjTiltedSpacalDetector::Construct_Fibers_SameLengthFiberPerTower::"
0683         << GetName() << " - constructed tower ID " << g_tower.id << " with "
0684         << fiber_count << " fibers. Average fiber length cut = "
0685         << std::accumulate(fiber_cut.begin(), fiber_cut.end(), 0.0) / fiber_cut.size() << " mm" << std::endl;
0686   }
0687 
0688   return fiber_count;
0689 }
0690 
0691 //! a block along z axis built with G4Trd that is slightly tapered in x dimension
0692 int PHG4FullProjTiltedSpacalDetector::Construct_Fibers(
0693     const PHG4FullProjTiltedSpacalDetector::SpacalGeom_t::geom_tower& g_tower,
0694     G4LogicalVolume* LV_tower)
0695 {
0696   G4Vector3D v_zshift = G4Vector3D(tan(g_tower.pTheta) * cos(g_tower.pPhi),
0697                                    tan(g_tower.pTheta) * sin(g_tower.pPhi), 1) *
0698                         g_tower.pDz;
0699   int fiber_cnt = 0;
0700   for (int ix = 0; ix < g_tower.NFiberX; ix++)
0701   {
0702     const double weighted_ix = static_cast<double>(ix) / (g_tower.NFiberX - 1.);
0703 
0704     const double weighted_pDx1 = (g_tower.pDx1 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0705     const double weighted_pDx2 = (g_tower.pDx2 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0706 
0707     const double weighted_pDx3 = (g_tower.pDx3 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0708     const double weighted_pDx4 = (g_tower.pDx4 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0709 
0710     for (int iy = 0; iy < g_tower.NFiberY; iy++)
0711     {
0712       if ((ix + iy) % 2 == 1)
0713       {
0714         continue;  // make a triangle pattern
0715       }
0716       const int fiber_ID = g_tower.compose_fiber_id(ix, iy);
0717 
0718       const double weighted_iy = static_cast<double>(iy) / (g_tower.NFiberY - 1.);
0719 
0720       const double weighted_pDy1 = (g_tower.pDy1 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
0721       const double weighted_pDy2 = (g_tower.pDy2 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
0722 
0723       const double weighted_pDx12 = weighted_pDx1 * (1 - weighted_iy) + weighted_pDx2 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp1);
0724       const double weighted_pDx34 = weighted_pDx3 * (1 - weighted_iy) + weighted_pDx4 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp2);
0725 
0726       G4Vector3D v1 = G4Vector3D(weighted_pDx12, weighted_pDy1, 0) - v_zshift;
0727       G4Vector3D v2 = G4Vector3D(weighted_pDx34, weighted_pDy2, 0) + v_zshift;
0728 
0729       G4Vector3D vector_fiber = (v2 - v1);
0730       vector_fiber *= (vector_fiber.mag() - get_geom_v3()->get_fiber_outer_r()) / vector_fiber.mag();  // shrink by fiber boundary protection
0731       G4Vector3D center_fiber = (v2 + v1) / 2;
0732 
0733       // convert to Geant4 units
0734       vector_fiber *= cm;
0735       center_fiber *= cm;
0736 
0737       const G4double fiber_length = vector_fiber.mag();
0738 
0739       std::stringstream ss;
0740       ss << std::string("_Tower") << g_tower.id;
0741       ss << "_x" << ix;
0742       ss << "_y" << iy;
0743       G4LogicalVolume* fiber_logic = Construct_Fiber(fiber_length,
0744                                                      ss.str());
0745 
0746       if (get_geom_v3()->get_construction_verbose() >= 3)
0747       {
0748         std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_Fibers::" << GetName()
0749                   << " - constructed fiber " << fiber_ID << ss.str()  //
0750                   << ", Length = " << fiber_length << "mm, "          //
0751                   << "x = " << center_fiber.x() << "mm, "             //
0752                   << "y = " << center_fiber.y() << "mm, "             //
0753                   << "z = " << center_fiber.z() << "mm, "             //
0754                   << "vx = " << vector_fiber.x() << "mm, "            //
0755                   << "vy = " << vector_fiber.y() << "mm, "            //
0756                   << "vz = " << vector_fiber.z() << "mm, "            //
0757                   << std::endl;
0758       }
0759 
0760       const G4double rotation_angle = G4Vector3D(0, 0, 1).angle(
0761           vector_fiber);
0762       const G4Vector3D rotation_axis =
0763           rotation_angle == 0 ? G4Vector3D(1, 0, 0) : G4Vector3D(0, 0, 1).cross(vector_fiber);
0764 
0765       G4Transform3D fiber_place(
0766           G4Translate3D(center_fiber.x(), center_fiber.y(),
0767                         center_fiber.z()) *
0768           G4Rotate3D(rotation_angle, rotation_axis));
0769 
0770       std::stringstream name;
0771       name << GetName() + std::string("_Tower") << g_tower.id << "_fiber"
0772            << ss.str();
0773 
0774       const bool overlapcheck_fiber = OverlapCheck() and (get_geom_v3()->get_construction_verbose() >= 3);
0775       G4PVPlacement* fiber_physi = new G4PVPlacement(fiber_place,
0776                                                      fiber_logic, G4String(name.str()), LV_tower, false,
0777                                                      fiber_ID, overlapcheck_fiber);
0778       fiber_vol[fiber_physi] = fiber_ID;
0779       assert(gdml_config);
0780       gdml_config->exclude_physical_vol(fiber_physi);
0781 
0782       ++fiber_cnt;
0783     }
0784   }
0785 
0786   if (get_geom_v3()->get_construction_verbose() >= 3)
0787   {
0788     std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_Fibers::" << GetName()
0789               << " - constructed tower ID " << g_tower.id << " with " << fiber_cnt
0790               << " fibers" << std::endl;
0791   }
0792 
0793   return fiber_cnt;
0794 }
0795 
0796 //! a block along z axis built with G4Trd that is slightly tapered in x dimension
0797 G4LogicalVolume*
0798 PHG4FullProjTiltedSpacalDetector::Construct_Tower(
0799     const PHG4FullProjTiltedSpacalDetector::SpacalGeom_t::geom_tower& g_tower)
0800 {
0801   std::stringstream sout;
0802   sout << "_" << g_tower.id;
0803   const G4String sTowerID(sout.str());
0804 
0805   // Processed PostionSeeds 1 from 1 1
0806 
0807   G4Trap* block_solid = new G4Trap(
0808       /*const G4String& pName*/ G4String(GetName()) + sTowerID,
0809       g_tower.pDz * cm,                                         // G4double pDz,
0810       g_tower.pTheta * rad, g_tower.pPhi * rad,                 // G4double pTheta, G4double pPhi,
0811       g_tower.pDy1 * cm, g_tower.pDx1 * cm, g_tower.pDx2 * cm,  // G4double pDy1, G4double pDx1, G4double pDx2,
0812       g_tower.pAlp1 * rad,                                      // G4double pAlp1,
0813       g_tower.pDy2 * cm, g_tower.pDx3 * cm, g_tower.pDx4 * cm,  // G4double pDy2, G4double pDx3, G4double pDx4,
0814       g_tower.pAlp2 * rad                                       // G4double pAlp2 //
0815   );
0816 
0817   G4Material* cylinder_mat = GetDetectorMaterial(get_geom_v3()->get_absorber_mat());
0818   assert(cylinder_mat);
0819 
0820   G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, cylinder_mat,
0821                                                      G4String(G4String(GetName()) + std::string("_Tower") + sTowerID), nullptr, nullptr,
0822                                                      nullptr);
0823 
0824   GetDisplayAction()->AddVolume(block_logic, "Block");
0825 
0826   // construct fibers
0827 
0828   if (get_geom_v3()->get_config() == SpacalGeom_t::kFullProjective_2DTaper_Tilted)
0829   {
0830     int fiber_count = Construct_Fibers(g_tower, block_logic);
0831 
0832     if (get_geom_v3()->get_construction_verbose() >= 2)
0833     {
0834       std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_Tower::" << GetName()
0835                 << " - constructed tower ID " << g_tower.id << " with "
0836                 << fiber_count << " fibers using Construct_Fibers" << std::endl;
0837     }
0838   }
0839   else if (get_geom_v3()->get_config() == SpacalGeom_t::kFullProjective_2DTaper_Tilted_SameLengthFiberPerTower)
0840   {
0841     int fiber_count = Construct_Fibers_SameLengthFiberPerTower(g_tower,
0842                                                                block_logic);
0843 
0844     if (get_geom_v3()->get_construction_verbose() >= 2)
0845     {
0846       std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_Tower::" << GetName()
0847                 << " - constructed tower ID " << g_tower.id << " with "
0848                 << fiber_count
0849                 << " fibers using Construct_Fibers_SameLengthFiberPerTower."
0850                 << "V = " << block_solid->GetCubicVolume() / (cm3) << "cm3, "
0851                 << "m = " << block_logic->GetMass() / gram << "gram, "
0852                 << "Density = " << (block_logic->GetMass() / gram) / (block_solid->GetCubicVolume() / cm3) << "g/cm3"
0853                 << std::endl;
0854     }
0855   }
0856   else
0857   {
0858     G4ExceptionDescription message;
0859     message << "can not recognize configuration type " << get_geom_v3()->get_config();
0860 
0861     G4Exception("PHG4FullProjTiltedSpacalDetector::Construct_Tower", "Wrong",
0862                 FatalException, message, "");
0863   }
0864 
0865   return block_logic;
0866 }
0867 
0868 G4LogicalVolume*
0869 PHG4FullProjTiltedSpacalDetector::Construct_LightGuide(
0870     const PHG4FullProjTiltedSpacalDetector::SpacalGeom_t::geom_tower& g_tower,
0871     const int index_x, const int index_y)
0872 {
0873   assert(_geom);
0874   std::stringstream sout;
0875   sout << "_Lightguide_" << g_tower.id << "_" << index_x << "_" << index_y;
0876   const G4String sTowerID(sout.str());
0877 
0878   assert(g_tower.LightguideHeight > 0);
0879 
0880   // light guide parameters in PHENIX units
0881   const double weight_x1 = 1 - (double) index_y / g_tower.NSubtowerY;
0882   const double weight_x2 = 1 - (double) (index_y + 1) / g_tower.NSubtowerY;
0883   const double weight_xcenter = 1 - (double) (index_y + 0.5) / g_tower.NSubtowerY;
0884 
0885   assert(weight_x1 >= 0 and weight_x1 <= 1);
0886   assert(weight_x2 >= 0 and weight_x2 <= 1);
0887   assert(weight_xcenter >= 0 and weight_xcenter <= 1);
0888 
0889   const double lg_pDx1 = (g_tower.pDx1 * weight_x1  //
0890                           + g_tower.pDx2 * (1 - weight_x1)) /
0891                          g_tower.NSubtowerX;
0892   const double lg_pDx2 = (g_tower.pDx1 * weight_x2  //
0893                           + g_tower.pDx2 * (1 - weight_x2)) /
0894                          g_tower.NSubtowerX;
0895   const double lg_pDy1 = g_tower.pDy1 / g_tower.NSubtowerY;
0896   const double lg_Alp1 = atan(
0897       (g_tower.pDx2 - g_tower.pDx1) * (-g_tower.NSubtowerX + 1. + 2 * index_x) / (double) (g_tower.NSubtowerX) / (2. * g_tower.pDy1) + tan(g_tower.pAlp1));
0898 
0899   const double shift_xcenter = (g_tower.pDx1 * weight_xcenter           //
0900                                 + g_tower.pDx2 * (1 - weight_xcenter))  //
0901                                *                                        //
0902                                (-g_tower.NSubtowerX + 1. + 2 * index_x) / (double) (g_tower.NSubtowerX);
0903   const double shift_ycenter = g_tower.pDy1  //
0904                                *             //
0905                                (-g_tower.NSubtowerY + 1. + 2 * index_y) / (double) (g_tower.NSubtowerY);
0906 
0907   G4VSolid* block_solid = new G4Trap(
0908       /*const G4String& pName*/ G4String(GetName()) + sTowerID,
0909       0.5 * g_tower.LightguideHeight * cm,  // G4double pDz,
0910       0 * rad, 0 * rad,                     // G4double pTheta, G4double pPhi,
0911       g_tower.LightguideTaperRatio * lg_pDy1 * cm,
0912       g_tower.LightguideTaperRatio * lg_pDx1 * cm,
0913       g_tower.LightguideTaperRatio * lg_pDx2 * cm,  // G4double pDy1, G4double pDx1, G4double pDx2,
0914       lg_Alp1 * rad,                                // G4double pAlp1,
0915       lg_pDy1 * cm, lg_pDx1 * cm, lg_pDx2 * cm,     // G4double pDy2, G4double pDx3, G4double pDx4,
0916       lg_Alp1 * rad                                 // G4double pAlp2 //
0917   );
0918 
0919   block_solid = new G4DisplacedSolid(G4String(GetName() + "_displaced"),
0920                                      block_solid, nullptr,                                     //
0921                                      G4ThreeVector(                                            //
0922                                          tan(g_tower.pTheta * rad) * cos(g_tower.pPhi * rad),  //
0923                                          tan(g_tower.pTheta * rad) * sin(g_tower.pPhi * rad),  //
0924                                          1) *                                                  // G4ThreeVector
0925                                              -(g_tower.pDz) *
0926                                              cm                                                       //
0927                                          + G4ThreeVector(shift_xcenter * cm, shift_ycenter * cm, 0)   // shit in subtower direction
0928                                          + G4ThreeVector(0, 0, -0.5 * g_tower.LightguideHeight * cm)  // shift in the light guide height
0929   );
0930 
0931   G4Material* cylinder_mat = GetDetectorMaterial(g_tower.LightguideMaterial);
0932   assert(cylinder_mat);
0933 
0934   G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, cylinder_mat,
0935                                                      G4String(G4String(GetName()) + std::string("_Tower") + sTowerID), nullptr, nullptr,
0936                                                      nullptr);
0937 
0938   GetDisplayAction()->AddMaterial("LightGuide", g_tower.LightguideMaterial);
0939   GetDisplayAction()->AddVolume(block_logic, "LightGuide");
0940 
0941   return block_logic;
0942 }
0943 
0944 void PHG4FullProjTiltedSpacalDetector::Print(const std::string& /*what*/) const
0945 {
0946   std::cout << "PHG4FullProjTiltedSpacalDetector::Print::" << GetName()
0947             << " - Print Geometry:" << std::endl;
0948   get_geom_v3()->Print();
0949 
0950   return;
0951 }
[ 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!