Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:20:45 

0001 #include "PHG4SpacalPrototype4Detector.h"
0002 
0003 #include <g4detectors/PHG4CylinderGeomContainer.h>
0004 #include <g4detectors/PHG4CylinderGeom_Spacalv1.h>  // for PHG4CylinderGeom_Spacalv1:...
0005 #include <g4detectors/PHG4CylinderGeom_Spacalv3.h>  // for PHG4CylinderGeom_...
0006 
0007 #include <phparameter/PHParameters.h>
0008 
0009 #include <g4main/PHG4Detector.h>  // for PHG4Detector
0010 #include <g4main/PHG4Utils.h>
0011 
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/PHIODataNode.h>
0014 #include <phool/PHNode.h>          // for PHNode
0015 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0016 #include <phool/PHObject.h>        // for PHObject
0017 #include <phool/getClass.h>
0018 
0019 #include <Geant4/G4Box.hh>
0020 #include <Geant4/G4Colour.hh>  // for G4Colour
0021 #include <Geant4/G4DisplacedSolid.hh>
0022 #include <Geant4/G4LogicalVolume.hh>
0023 #include <Geant4/G4Material.hh>
0024 #include <Geant4/G4PVPlacement.hh>
0025 #include <Geant4/G4PhysicalConstants.hh>
0026 #include <Geant4/G4String.hh>  // for G4String
0027 #include <Geant4/G4SubtractionSolid.hh>
0028 #include <Geant4/G4SystemOfUnits.hh>
0029 #include <Geant4/G4ThreeVector.hh>  // for G4ThreeVector
0030 #include <Geant4/G4Transform3D.hh>  // for G4Transform3D, G4Translate3D
0031 #include <Geant4/G4Trap.hh>
0032 #include <Geant4/G4Tubs.hh>
0033 #include <Geant4/G4Types.hh>  // for G4double
0034 #include <Geant4/G4UserLimits.hh>
0035 #include <Geant4/G4Vector3D.hh>
0036 #include <Geant4/G4VisAttributes.hh>
0037 
0038 #include <boost/foreach.hpp>
0039 
0040 #include <algorithm>
0041 #include <cassert>
0042 #include <cmath>
0043 #include <cstdlib>   // for exit
0044 #include <iostream>  // for operator<<, basic_ostream
0045 #include <limits>    // for numeric_limits
0046 #include <memory>    // for allocator_traits<...
0047 #include <numeric>   // std::accumulate
0048 #include <sstream>
0049 #include <string>  // std::string, std::to_string
0050 #include <vector>  // for vector
0051 
0052 class PHG4CylinderGeom;
0053 
0054 using namespace std;
0055 
0056 //_______________________________________________________________
0057 //note this inactive thickness is ~1.5% of a radiation length
0058 PHG4SpacalPrototype4Detector::PHG4SpacalPrototype4Detector(PHG4Subsystem* subsys, PHCompositeNode* Node, PHParameters* parameters, const std::string& dnam)
0059   : PHG4Detector(subsys, Node, dnam)
0060   , construction_params(parameters)
0061   , cylinder_solid(nullptr)
0062   , cylinder_logic(nullptr)
0063   , cylinder_physi(nullptr)
0064   ,  //
0065   active(0)
0066   , absorberactive(0)
0067   ,  //
0068   step_limits(nullptr)
0069   , clading_step_limits(nullptr)
0070   , fiber_core_step_limits(nullptr)
0071   ,  //
0072   _geom(nullptr)
0073 {
0074   SetActive(construction_params->get_int_param("active"));
0075   SetAbsorberActive(construction_params->get_int_param("absorberactive"));
0076   SuperDetector(dnam);
0077 }
0078 
0079 PHG4SpacalPrototype4Detector::~PHG4SpacalPrototype4Detector(void)
0080 {
0081   // deleting nullptr pointers is allowed (results in NOOP)
0082   // so checking for not null before deleting is not needed
0083   delete step_limits;
0084   delete clading_step_limits;
0085   delete fiber_core_step_limits;
0086   delete _geom;
0087 }
0088 
0089 //_______________________________________________________________
0090 int PHG4SpacalPrototype4Detector::IsInCylinderActive(
0091     const G4VPhysicalVolume* volume)
0092 {
0093   //  cout << "checking detector" << endl;
0094   if (active && fiber_core_vol.find(volume) != fiber_core_vol.end())
0095   {
0096     //      return fiber_core_vol.find(volume)->second;
0097     return FIBER_CORE;
0098   }
0099   if (absorberactive)
0100   {
0101     if (fiber_vol.find(volume) != fiber_vol.end())
0102       return FIBER_CLADING;
0103 
0104     if (block_vol.find(volume) != block_vol.end())
0105       return ABSORBER;
0106 
0107     if (calo_vol.find(volume) != calo_vol.end())
0108       return SUPPORT;
0109   }
0110   return INACTIVE;
0111 }
0112 
0113 //_______________________________________________________________
0114 void PHG4SpacalPrototype4Detector::ConstructMe(G4LogicalVolume* logicWorld)
0115 {
0116   PHNodeIterator iter(topNode());
0117   PHCompositeNode* parNode = dynamic_cast<PHCompositeNode*>(iter.findFirst(
0118       "PHCompositeNode", "RUN"));
0119   assert(parNode);
0120 
0121   if (!_geom)
0122   {
0123     _geom = new SpacalGeom_t();
0124   }
0125   assert(_geom);
0126 
0127   _geom->set_config(
0128       SpacalGeom_t::kFullProjective_2DTaper_Tilted_SameLengthFiberPerTower);
0129   //  _geom->load_demo_sector_tower_map4();
0130   //    _geom->set_construction_verbose(2);
0131   _geom->ImportParameters(*construction_params);
0132 
0133   _geom->subtower_consistency_check();
0134 
0135   //  step_limits = new G4UserLimits(_geom->get_calo_step_size() * cm);
0136   //
0137   //  clading_step_limits = new G4UserLimits(
0138   //      _geom->get_fiber_clading_step_size() * cm);
0139   step_limits = nullptr;
0140   clading_step_limits = nullptr;
0141 
0142   fiber_core_step_limits = new G4UserLimits(
0143       _geom->get_fiber_core_step_size() * cm);
0144 
0145   const G4double z_rotation = construction_params->get_double_param(
0146                                   "z_rotation_degree") *
0147                               degree;
0148   const G4double enclosure_x = construction_params->get_double_param(
0149                                    "enclosure_x") *
0150                                cm;
0151   const G4double enclosure_x_shift = construction_params->get_double_param(
0152                                          "enclosure_x_shift") *
0153                                      cm;
0154   const G4double enclosure_y = construction_params->get_double_param(
0155                                    "enclosure_y") *
0156                                cm;
0157   const G4double enclosure_z = construction_params->get_double_param(
0158                                    "enclosure_z") *
0159                                cm;
0160 
0161   const G4double box_x_shift = (_geom->get_radius() + 0.5 * _geom->get_thickness()) * cm + enclosure_x_shift;
0162   const G4double box_z_shift = 0.5 * (_geom->get_zmin() + _geom->get_zmax()) * cm;
0163 
0164   //  G4Tubs* _cylinder_solid = new G4Tubs(G4String(GetName()),
0165   //      _geom->get_radius() * cm, _geom->get_max_radius() * cm,
0166   //      _geom->get_length() * cm / 2.0, 0, twopi);
0167   //  G4VSolid* cylinder_solid = new G4Box(G4String(GetName()),
0168   //      _geom->get_thickness() * 1.1 * 0.5 * cm, //
0169   //      twopi / _geom->get_azimuthal_n_sec() * _geom->get_sector_map().size()
0170   //          * 0.5 * _geom->get_max_radius() * cm, //
0171   //      _geom->get_length() * cm / 2.0);
0172   G4VSolid* cylinder_solid = new G4Box(G4String(GetName()),
0173                                        enclosure_x * 0.5,  //
0174                                        enclosure_y * 0.5,  //
0175                                        enclosure_z * 0.5);
0176 
0177   cylinder_solid = new G4DisplacedSolid(G4String(GetName() + "_displaced"),
0178                                         cylinder_solid, 0,  //
0179                                         G4ThreeVector(box_x_shift, 0, box_z_shift));
0180 
0181   G4Material* cylinder_mat = G4Material::GetMaterial("G4_AIR");
0182   assert(cylinder_mat);
0183 
0184   cylinder_logic = new G4LogicalVolume(cylinder_solid, cylinder_mat,
0185                                        G4String(GetName()), 0, 0, 0);
0186   G4VisAttributes* VisAtt = new G4VisAttributes();
0187   VisAtt->SetColour(G4Colour::White());
0188   VisAtt->SetVisibility(true);
0189   VisAtt->SetForceSolid(false);
0190   VisAtt->SetForceWireframe(true);
0191   cylinder_logic->SetVisAttributes(VisAtt);
0192 
0193   G4Transform3D cylinder_place(
0194       G4Translate3D(_geom->get_xpos() * cm, _geom->get_ypos() * cm,
0195                     _geom->get_zpos() * cm)  //
0196       * G4RotateZ3D(z_rotation)              //
0197       * G4Translate3D(
0198             -(_geom->get_radius() + 0.5 * _geom->get_thickness()) * cm, 0,
0199             0));
0200 
0201   cylinder_physi = new G4PVPlacement(cylinder_place, cylinder_logic,
0202                                      G4String(GetName()), logicWorld, false, 0, OverlapCheck());
0203 
0204   // install sectors
0205   std::pair<G4LogicalVolume*, G4Transform3D> psec = Construct_AzimuthalSeg();
0206   G4LogicalVolume* sec_logic = psec.first;
0207   const G4Transform3D& sec_trans = psec.second;
0208   BOOST_FOREACH (const SpacalGeom_t::sector_map_t::value_type& val, _geom->get_sector_map())
0209   {
0210     const int sec = val.first;
0211     const double rot = val.second;
0212 
0213     G4Transform3D sec_place = G4RotateZ3D(rot) * sec_trans;
0214 
0215     ostringstream name;
0216     name << GetName() << "_sec" << sec;
0217 
0218     G4PVPlacement* calo_phys = new G4PVPlacement(sec_place, sec_logic,
0219                                                  G4String(name.str()), cylinder_logic, false, sec,
0220                                                  OverlapCheck());
0221     calo_vol[calo_phys] = sec;
0222   }
0223   _geom->set_nscint(_geom->get_nscint() * _geom->get_sector_map().size());
0224 
0225   // install electronics
0226   const G4double electronics_thickness = construction_params->get_double_param(
0227                                              "electronics_thickness") *
0228                                          cm;
0229   const string electronics_material = construction_params->get_string_param(
0230       "electronics_material");
0231 
0232   G4VSolid* electronics_solid = new G4Box(G4String(GetName()),
0233                                           electronics_thickness / 2.0,  //
0234                                           sin(
0235                                               twopi / _geom->get_azimuthal_n_sec() * _geom->get_sector_map().size() / 2) *
0236                                               (_geom->get_radius() - _geom->get_max_lightguide_height()) * cm,  //
0237                                           _geom->get_length() * cm / 2.0);
0238 
0239   electronics_solid = new G4DisplacedSolid(G4String(GetName() + "_displaced"),
0240                                            electronics_solid,
0241                                            0,  //
0242                                            G4ThreeVector(
0243                                                cos(
0244                                                    twopi / _geom->get_azimuthal_n_sec() * _geom->get_sector_map().size() / 2) *
0245                                                        (_geom->get_radius() - _geom->get_max_lightguide_height()) * cm -
0246                                                    electronics_thickness,
0247                                                0, box_z_shift));
0248 
0249   G4Material* electronics_mat = G4Material::GetMaterial(electronics_material);
0250   assert(electronics_mat);
0251 
0252   G4LogicalVolume* electronics_logic = new G4LogicalVolume(electronics_solid,
0253                                                            electronics_mat, G4String(GetName()) + "_Electronics", 0, 0, 0);
0254   VisAtt = new G4VisAttributes();
0255   PHG4Utils::SetColour(VisAtt, electronics_material);
0256   VisAtt->SetVisibility(true);
0257   VisAtt->SetForceSolid(not _geom->is_virualize_fiber());
0258   electronics_logic->SetVisAttributes(VisAtt);
0259 
0260   new G4PVPlacement(G4Transform3D::Identity, electronics_logic,
0261                     G4String(GetName()) + "_Electronics", cylinder_logic, false, 0,
0262                     OverlapCheck());
0263 
0264   // install the enclosure box
0265   const G4double enclosure_thickness = construction_params->get_double_param(
0266                                            "enclosure_thickness") *
0267                                        cm;
0268   const string enclosure_material = construction_params->get_string_param(
0269       "enclosure_material");
0270 
0271   G4VSolid* outer_enclosur_solid = new G4Box(
0272       G4String(GetName()) + "_outer_enclosur_solid", enclosure_x * 0.5,  //
0273       enclosure_y * 0.5,                                                 //
0274       enclosure_z * 0.5);
0275   G4VSolid* inner_enclosur_solid = new G4Box(
0276       G4String(GetName()) + "_inner_enclosur_solid",
0277       enclosure_x * 0.5 - enclosure_thickness,  //
0278       enclosure_y * 0.5 - enclosure_thickness,  //
0279       enclosure_z * 0.5 - enclosure_thickness);
0280   G4VSolid* enclosure_solid = new G4SubtractionSolid(
0281       G4String(GetName()) + "_enclosure_solid",  //
0282       outer_enclosur_solid, inner_enclosur_solid);
0283   enclosure_solid = new G4DisplacedSolid(
0284       G4String(GetName() + "_enclosure_solid_displaced"), enclosure_solid, 0,  //
0285       G4ThreeVector(box_x_shift, 0, box_z_shift));
0286 
0287   G4Material* enclosure_mat = G4Material::GetMaterial(enclosure_material);
0288   assert(enclosure_mat);
0289 
0290   G4LogicalVolume* enclosure_logic = new G4LogicalVolume(enclosure_solid,
0291                                                          enclosure_mat, G4String(GetName()) + "_enclosure", 0, 0, 0);
0292   VisAtt = new G4VisAttributes();
0293   PHG4Utils::SetColour(VisAtt, enclosure_material);
0294   VisAtt->SetVisibility(true);
0295   VisAtt->SetForceSolid(not _geom->is_virualize_fiber());
0296   enclosure_logic->SetVisAttributes(VisAtt);
0297 
0298   new G4PVPlacement(G4Transform3D::Identity, enclosure_logic,
0299                     G4String(GetName()) + "_enclosure", cylinder_logic, false, 0,
0300                     OverlapCheck());
0301 
0302   if (active)
0303   {
0304     ostringstream geonode;
0305     if (superdetector != "NONE")
0306     {
0307       geonode << "CYLINDERGEOM_" << superdetector;
0308     }
0309     else
0310     {
0311       geonode << "CYLINDERGEOM_" << detector_type;
0312     }
0313     PHG4CylinderGeomContainer* geo = findNode::getClass<
0314         PHG4CylinderGeomContainer>(topNode(), geonode.str());
0315     if (!geo)
0316     {
0317       geo = new PHG4CylinderGeomContainer();
0318       PHNodeIterator iter(topNode());
0319       PHCompositeNode* runNode =
0320           dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode",
0321                                                         "RUN"));
0322       PHIODataNode<PHObject>* newNode = new PHIODataNode<PHObject>(geo,
0323                                                                    geonode.str(), "PHObject");
0324       runNode->addNode(newNode);
0325     }
0326     // here in the detector class we have internal units, convert to cm
0327     // before putting into the geom object
0328     PHG4CylinderGeom* mygeom = new SpacalGeom_t(*_geom);
0329     geo->AddLayerGeom(0, mygeom);
0330     if (_geom->get_construction_verbose() >= 1)
0331     {
0332       cout << "PHG4SpacalPrototype4Detector::Construct::" << GetName()
0333            << " - Print Layer Geometry:" << endl;
0334       geo->identify();
0335     }
0336   }
0337 
0338   if (absorberactive)
0339   {
0340     ostringstream geonode;
0341     if (superdetector != "NONE")
0342     {
0343       geonode << "CYLINDERGEOM_ABSORBER_" << superdetector;
0344     }
0345     else
0346     {
0347       geonode << "CYLINDERGEOM_ABSORBER_" << detector_type << "_" << 0;
0348     }
0349     PHG4CylinderGeomContainer* geo = findNode::getClass<PHG4CylinderGeomContainer>(topNode(), geonode.str());
0350     if (!geo)
0351     {
0352       geo = new PHG4CylinderGeomContainer();
0353       PHNodeIterator iter(topNode());
0354       PHCompositeNode* runNode =
0355           dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode",
0356                                                         "RUN"));
0357       PHIODataNode<PHObject>* newNode = new PHIODataNode<PHObject>(geo, geonode.str(), "PHObject");
0358       runNode->addNode(newNode);
0359     }
0360     // here in the detector class we have internal units, convert to cm
0361     // before putting into the geom object
0362     PHG4CylinderGeom* mygeom = new SpacalGeom_t(*_geom);
0363     geo->AddLayerGeom(0, mygeom);
0364     if (_geom->get_construction_verbose() >= 1)
0365     {
0366       cout << "PHG4SpacalPrototype4Detector::Construct::" << GetName()
0367            << " - Print Absorber Layer Geometry:" << endl;
0368       geo->identify();
0369     }
0370   }
0371 
0372   if (_geom->get_construction_verbose() >= 1)
0373   {
0374     cout << "PHG4SpacalPrototype4Detector::Construct::" << GetName()
0375          << " - Completed. Print G4 Geometry:" << endl;
0376     Print();
0377     cout << "box_x_shift = " << box_x_shift << endl;
0378     cout << "box_z_shift = " << box_z_shift << endl;
0379   }
0380 }
0381 
0382 std::pair<G4LogicalVolume*, G4Transform3D>
0383 PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg()
0384 {
0385   assert(_geom);
0386   assert(_geom->get_azimuthal_n_sec() > 4);
0387 
0388   // basic tilt geometry
0389   const G4double half_chord_backend =
0390       _geom->get_max_radius() * cm * tan(pi / _geom->get_azimuthal_n_sec())  //
0391       + fabs(_geom->get_thickness() * cm * 0.5 * tan(_geom->get_azimuthal_tilt()));
0392   const G4double reduced_outer_radius = sqrt(pow(_geom->get_max_radius() * cm, 2) - half_chord_backend * half_chord_backend);
0393   const G4double enclosure_depth = reduced_outer_radius - _geom->get_radius() * cm;
0394   const G4double enclosure_center = 0.5 * (reduced_outer_radius + _geom->get_radius() * cm);
0395   const G4double enclosure_half_height_half_width = enclosure_center * tan(pi / _geom->get_azimuthal_n_sec());
0396 
0397   const G4double width_adj1 = tan(_geom->get_azimuthal_tilt() - pi / _geom->get_azimuthal_n_sec()) * enclosure_depth * 0.5;
0398   const G4double width_adj2 = tan(_geom->get_azimuthal_tilt() + pi / _geom->get_azimuthal_n_sec()) * enclosure_depth * 0.5;
0399 
0400   const G4double center_adj = (width_adj1 + width_adj2) * 0.5;
0401   const G4double center_tilt_angle = atan2(center_adj, enclosure_depth * 0.5);
0402   const G4double inner_half_width = enclosure_half_height_half_width + 0.5 * (width_adj1 - width_adj2);
0403   const G4double outter_half_width = enclosure_half_height_half_width + 0.5 * (-width_adj1 + width_adj2);
0404 
0405   // enclosure walls
0406   const G4double edge1_tilt_angle = atan2(width_adj1, enclosure_depth * 0.5);
0407   const G4double edge2_tilt_angle = atan2(width_adj2, enclosure_depth * 0.5);
0408   //  const G4double edge1_half_depth = sqrt(width_adj1 * width_adj1 + enclosure_depth * enclosure_depth * .25);
0409   //  const G4double edge2_half_depth = sqrt(width_adj2 * width_adj2 + enclosure_depth * enclosure_depth * .25);
0410 
0411   // projective center
0412   const G4double half_projection_ratio = 0.5 * (-width_adj1 + width_adj2) / enclosure_half_height_half_width;
0413   const G4double projection_center_y = enclosure_center - ((enclosure_depth * 0.5) / half_projection_ratio);
0414   const G4double projection_center_x = center_adj / half_projection_ratio;
0415 
0416   // blocks azimuthal segmentation
0417   const int phi_bin_in_sec = _geom->get_max_phi_bin_in_sec();
0418   assert(phi_bin_in_sec >= 1);
0419   const G4double block_azimuth_angle = (edge2_tilt_angle - edge1_tilt_angle) / phi_bin_in_sec;
0420   assert(block_azimuth_angle > 0);
0421   if (!(fabs(block_azimuth_angle - M_PI * 2 / _geom->get_azimuthal_n_sec() / phi_bin_in_sec) < M_PI * numeric_limits<G4double>::epsilon()))
0422   {
0423     cout << "angle/nsec out of range: " << M_PI * numeric_limits<G4double>::epsilon() << endl;
0424     exit(1);
0425   }
0426   const G4double block_edge1_half_width = enclosure_half_height_half_width - (_geom->get_sidewall_thickness() * cm + _geom->get_sidewall_outer_torr() * cm + 2.0 * _geom->get_assembly_spacing() * cm) / cos(edge1_tilt_angle);
0427   const G4double block_edge2_half_width = enclosure_half_height_half_width - (_geom->get_sidewall_thickness() * cm + _geom->get_sidewall_outer_torr() * cm + 2.0 * _geom->get_assembly_spacing() * cm) / cos(edge2_tilt_angle);
0428   G4double block_width_ratio = 0;
0429   for (int s = 0; s < phi_bin_in_sec; ++s)
0430   {
0431     block_width_ratio += 1 / cos(block_azimuth_angle * (0.5 + s) + edge1_tilt_angle);
0432   }
0433   const G4double block_half_height_width = (block_edge1_half_width + block_edge2_half_width) / block_width_ratio;
0434   assert(block_half_height_width > 0);
0435 
0436   // write out the azimuthal block geometry
0437   // block azimuth geometry records
0438   struct block_azimuth_geom
0439   {
0440     G4double angle;
0441     G4double projection_center_y;
0442     G4double projection_center_x;
0443     G4double projection_length;
0444   };
0445   vector<block_azimuth_geom> block_azimuth_geoms(phi_bin_in_sec,
0446                                                  block_azimuth_geom{
0447                                                      numeric_limits<double>::signaling_NaN(),
0448                                                      numeric_limits<double>::signaling_NaN(),
0449                                                      numeric_limits<double>::signaling_NaN(),
0450                                                      numeric_limits<double>::signaling_NaN()});  // [phi-bin in sector] -> azimuth geometry
0451   G4double block_x_edge1 = block_edge1_half_width;
0452   for (int s = 0; s < phi_bin_in_sec; ++s)
0453   {
0454     block_azimuth_geom& geom = block_azimuth_geoms[s];
0455 
0456     geom.angle = block_azimuth_angle * (0.5 + s) + edge1_tilt_angle;
0457     const G4double block_x_size = block_half_height_width / cos(geom.angle);
0458     assert(block_x_size > 0);
0459     const G4double x_center = block_x_edge1 - 0.5 * block_x_size;
0460 
0461     // projection center per block
0462     geom.projection_length = block_half_height_width / 2. / tan(block_azimuth_angle / 2.);
0463     assert(geom.projection_length > 0);
0464     geom.projection_center_y = enclosure_center - geom.projection_length * cos(geom.angle);
0465     geom.projection_center_x = x_center + geom.projection_length * sin(geom.angle);
0466 
0467     // next step
0468     block_x_edge1 -= block_x_size;
0469   }
0470 
0471   //write out the azimuthal block divider's geometry
0472   struct block_divider_azimuth_geom
0473   {
0474     G4double angle;  //! rotation angle
0475     G4double projection_center_y;
0476     G4double projection_center_x;
0477     G4double thickness;            // thickness in the approximate azimuth direction
0478     G4double radial_displacement;  //! displacement along the width direction, which is the radial direction if tilt = 0
0479     G4double width;                //! wdith along the approximate radial direction
0480   };
0481   assert(phi_bin_in_sec >= 1);
0482   vector<block_divider_azimuth_geom> divider_azimuth_geoms(phi_bin_in_sec - 1,
0483                                                            block_divider_azimuth_geom{
0484                                                                numeric_limits<double>::signaling_NaN(),
0485                                                                numeric_limits<double>::signaling_NaN(),
0486                                                                numeric_limits<double>::signaling_NaN(),
0487                                                                numeric_limits<double>::signaling_NaN(),
0488                                                                numeric_limits<double>::signaling_NaN(),
0489                                                                numeric_limits<double>::signaling_NaN()});
0490 
0491   if (_geom->get_sidewall_thickness() > 0)
0492   {
0493     for (int s = 0; s < phi_bin_in_sec - 1; ++s)
0494     {
0495       block_divider_azimuth_geom& geom = divider_azimuth_geoms[s];
0496 
0497       geom.angle = 0.5 * (block_azimuth_geoms[s].angle + block_azimuth_geoms[s + 1].angle);
0498       geom.projection_center_y = 0.5 * (block_azimuth_geoms[s].projection_center_y + block_azimuth_geoms[s + 1].projection_center_y);
0499       geom.projection_center_x = 0.5 * (block_azimuth_geoms[s].projection_center_x + block_azimuth_geoms[s + 1].projection_center_x);
0500       geom.radial_displacement = 0.5 * (block_azimuth_geoms[s].projection_length + block_azimuth_geoms[s + 1].projection_length);
0501 
0502       geom.thickness = 2.0 * _geom->get_assembly_spacing() * cm * cos(block_azimuth_angle / 2.) - 2 * um;
0503       geom.width = _geom->get_divider_width() * cm;
0504     }
0505   }
0506 
0507   if (fabs(block_x_edge1 - (-block_edge2_half_width)) > _geom->get_assembly_spacing() * cm)
0508   {
0509     cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg - ERROR - " << endl
0510          << "\t block_x_edge1 = " << block_x_edge1 << endl
0511          << "\t block_edge2_half_width = " << block_edge2_half_width << endl
0512          << "\t fabs(block_x_edge1 - (-block_edge2_half_width)) = " << fabs(block_x_edge1 - (-block_edge2_half_width)) << endl
0513          << "\t _geom->get_assembly_spacing() * cm = " << _geom->get_assembly_spacing() * cm << endl;
0514   }
0515   if (!(fabs(block_x_edge1 - (-block_edge2_half_width)) < _geom->get_assembly_spacing() * cm))  // closure check
0516   {
0517     cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg - ERROR - closure check failed: " << fabs(block_x_edge1 - (-block_edge2_half_width)) << endl;
0518     exit(1);
0519   }
0520 
0521   if (Verbosity())
0522   {
0523     cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg - " << endl
0524          << "\t edge1_tilt_angle = " << edge1_tilt_angle << endl
0525          << "\t edge2_tilt_angle = " << edge2_tilt_angle << endl
0526          << "\t projection_center_y = " << projection_center_y << endl
0527          << "\t projection_center_x = " << projection_center_x << endl
0528          << "\t block_azimuth_angle = " << block_azimuth_angle << endl
0529          << "\t block_edge1_half_width = " << block_edge1_half_width << endl
0530          << "\t block_edge2_half_width = " << block_edge2_half_width << endl
0531          << "\t block_width_ratio = " << block_width_ratio << endl
0532          << "\t block_half_height_width = " << block_half_height_width << endl;
0533 
0534     for (int s = 0; s < phi_bin_in_sec; ++s)
0535     {
0536       cout << "\t block[" << s << "].angle = " << block_azimuth_geoms[s].angle << endl;
0537       cout << "\t block[" << s << "].projection_center_y = " << block_azimuth_geoms[s].projection_center_y << endl;
0538       cout << "\t block[" << s << "].projection_center_x = " << block_azimuth_geoms[s].projection_center_x << endl;
0539     }
0540     for (int s = 0; s < phi_bin_in_sec - 1; ++s)
0541     {
0542       cout << "\t divider[" << s << "].angle = " << divider_azimuth_geoms[s].angle << endl;
0543       cout << "\t divider[" << s << "].projection_center_x = " << divider_azimuth_geoms[s].projection_center_x << endl;
0544       cout << "\t divider[" << s << "].projection_center_y = " << divider_azimuth_geoms[s].projection_center_y << endl;
0545       cout << "\t divider[" << s << "].radial_displacement = " << divider_azimuth_geoms[s].radial_displacement << endl;
0546       cout << "\t divider[" << s << "].thickness = " << divider_azimuth_geoms[s].thickness << endl;
0547       cout << "\t divider[" << s << "].width = " << divider_azimuth_geoms[s].width << endl;
0548     }
0549   }
0550 
0551   assert(enclosure_depth > 10 * cm);
0552 
0553   G4VSolid* sec_solid = new G4Trap(
0554       /*const G4String& pName*/ G4String(GetName() + string("_sec_trap")),
0555       enclosure_depth * 0.5,                                                              // G4double pDz,
0556       center_tilt_angle, halfpi,                                                          // G4double pTheta, G4double pPhi,
0557       inner_half_width, _geom->get_length() * cm / 2.0, _geom->get_length() * cm / 2.0,   // G4double pDy1, G4double pDx1, G4double pDx2,
0558       0,                                                                                  // G4double pAlp1,
0559       outter_half_width, _geom->get_length() * cm / 2.0, _geom->get_length() * cm / 2.0,  // G4double pDy2, G4double pDx3, G4double pDx4,
0560       0                                                                                   // G4double pAlp2 //
0561   );
0562   const G4double box_z_shift = 0.5 * (_geom->get_zmin() + _geom->get_zmax()) * cm;
0563   G4Transform3D sec_solid_transform =
0564       G4TranslateZ3D(box_z_shift) * G4TranslateY3D(enclosure_center) * G4RotateY3D(halfpi) * G4RotateX3D(-halfpi);
0565   G4VSolid* sec_solid_place = new G4DisplacedSolid(
0566       G4String(GetName() + string("_sec")), sec_solid, sec_solid_transform);
0567 
0568   G4Material* cylinder_mat = G4Material::GetMaterial("G4_AIR");
0569   assert(cylinder_mat);
0570 
0571   G4LogicalVolume* sec_logic = new G4LogicalVolume(sec_solid_place, cylinder_mat,
0572                                                    G4String(G4String(GetName() + string("_sec"))), 0, 0, nullptr);
0573 
0574   G4VisAttributes* VisAtt = new G4VisAttributes();
0575   VisAtt->SetColor(.5, .9, .5, .1);
0576   VisAtt->SetVisibility(true);
0577   VisAtt->SetForceSolid(false);
0578   VisAtt->SetForceWireframe(true);
0579   sec_logic->SetVisAttributes(VisAtt);
0580 
0581   //  // construct walls
0582   //
0583   //  G4Material* wall_mat = G4Material::GetMaterial(_geom->get_sidewall_mat());
0584   //  assert(wall_mat);
0585   //
0586   //  if (_geom->get_sidewall_thickness() > 0)
0587   //  {
0588   //    // end walls
0589   //    if (_geom->get_construction_verbose() >= 1)
0590   //    {
0591   //      cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg::" << GetName()
0592   //           << " - construct end walls." << endl;
0593   //    }
0594   //    //        G4Tubs* wall_solid = new G4Tubs(G4String(GetName() + string("_EndWall")),
0595   //    //            _geom->get_radius() * cm + _geom->get_sidewall_outer_torr() * cm,
0596   //    //            _geom->get_max_radius() * cm - _geom->get_sidewall_outer_torr() * cm,
0597   //    //            _geom->get_sidewall_thickness() * cm / 2.0,
0598   //    //            halfpi - pi / _geom->get_azimuthal_n_sec(),
0599   //    //            twopi / _geom->get_azimuthal_n_sec());
0600   //    const G4double side_wall_half_thickness = _geom->get_sidewall_thickness() * cm / 2.0;
0601   //    G4VSolid* wall_solid = new G4Trap(G4String(GetName() + string("_EndWall_trap")),
0602   //                                      enclosure_depth * 0.5,                                                  // G4double pDz,
0603   //                                      center_tilt_angle, halfpi,                                              // G4double pTheta, G4double pPhi,
0604   //                                      inner_half_width, side_wall_half_thickness, side_wall_half_thickness,   // G4double pDy1, G4double pDx1, G4double pDx2,
0605   //                                      0,                                                                      // G4double pAlp1,
0606   //                                      outter_half_width, side_wall_half_thickness, side_wall_half_thickness,  // G4double pDy2, G4double pDx3, G4double pDx4,
0607   //                                      0                                                                       // G4double pAlp2 //
0608   //    );
0609   //    G4VSolid* wall_solid_place = new G4DisplacedSolid(
0610   //        G4String(GetName() + string("_EndWall")), wall_solid, sec_solid_transform);
0611   //
0612   //    G4LogicalVolume* wall_logic = new G4LogicalVolume(wall_solid_place, wall_mat,
0613   //                                                      G4String(G4String(GetName() + string("_EndWall"))), 0, 0,
0614   //                                                      nullptr);
0615   //    GetDisplayAction()->AddVolume(wall_logic, "Wall");
0616   //
0617   //    typedef map<int, double> z_locations_t;
0618   //    z_locations_t z_locations;
0619   //    z_locations[000] = _geom->get_sidewall_thickness() * cm / 2.0 + _geom->get_assembly_spacing() * cm;
0620   //    z_locations[001] = _geom->get_length() * cm / 2.0 - (_geom->get_sidewall_thickness() * cm / 2.0 + _geom->get_assembly_spacing() * cm);
0621   //    z_locations[100] = -(_geom->get_sidewall_thickness() * cm / 2.0 + _geom->get_assembly_spacing() * cm);
0622   //    z_locations[101] = -(_geom->get_length() * cm / 2.0 - (_geom->get_sidewall_thickness() * cm / 2.0 + _geom->get_assembly_spacing() * cm));
0623   //
0624   //    BOOST_FOREACH (z_locations_t::value_type& val, z_locations)
0625   //    {
0626   //      if (_geom->get_construction_verbose() >= 2)
0627   //        cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg::"
0628   //             << GetName() << " - constructed End Wall ID " << val.first
0629   //             << " @ Z = " << val.second << endl;
0630   //
0631   //      G4Transform3D wall_trans = G4TranslateZ3D(val.second);
0632   //
0633   //      G4PVPlacement* wall_phys = new G4PVPlacement(wall_trans, wall_logic,
0634   //                                                   G4String(GetName()) + G4String("_EndWall_") + to_string(val.first), sec_logic,
0635   //                                                   false, val.first, OverlapCheck());
0636   //
0637   //      calo_vol[wall_phys] = val.first;
0638   //      assert(gdml_config);
0639   //      gdml_config->exclude_physical_vol(wall_phys);
0640   //    }
0641   //  }
0642   //  //
0643   //  if (_geom->get_sidewall_thickness() > 0)
0644   //  {
0645   //    // side walls
0646   //    if (_geom->get_construction_verbose() >= 1)
0647   //    {
0648   //      cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg::" << GetName()
0649   //           << " - construct side walls." << endl;
0650   //    }
0651   //
0652   //    typedef map<int, pair<int, int> > sign_t;
0653   //    sign_t signs;
0654   //    signs[100] = make_pair(+1, +1);
0655   //    signs[101] = make_pair(+1, -1);
0656   //    signs[200] = make_pair(-1, +1);
0657   //    signs[201] = make_pair(-1, -1);
0658   //
0659   //    BOOST_FOREACH (sign_t::value_type& val, signs)
0660   //    {
0661   //      const int sign_z = val.second.first;
0662   //      const int sign_azimuth = val.second.second;
0663   //
0664   //      if (_geom->get_construction_verbose() >= 2)
0665   //        cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg::"
0666   //             << GetName() << " - constructed Side Wall ID " << val.first
0667   //             << " with"
0668   //             << " Shift X = "
0669   //             << sign_azimuth * (_geom->get_sidewall_thickness() * cm / 2.0 + _geom->get_sidewall_outer_torr() * cm)
0670   //             << " Rotation Z = "
0671   //             << sign_azimuth * pi / _geom->get_azimuthal_n_sec()
0672   //             << " Shift Z = " << sign_z * (_geom->get_length() * cm / 4)
0673   //             << endl;
0674   //
0675   //      const G4double azimuth_roate =
0676   //          sign_azimuth > 0 ? edge1_tilt_angle : edge2_tilt_angle;
0677   //      const G4double edge_half_depth = -_geom->get_sidewall_thickness() * cm - _geom->get_sidewall_outer_torr() * cm +
0678   //                                       (sign_azimuth > 0 ? edge1_half_depth : edge2_half_depth);
0679   //
0680   //      G4Box* wall_solid = new G4Box(G4String(GetName() + G4String("_SideWall_") + to_string(val.first)),
0681   //                                    _geom->get_sidewall_thickness() * cm / 2.0,
0682   //                                    edge_half_depth,
0683   //                                    (_geom->get_length() / 2. - 2 * (_geom->get_sidewall_thickness() + 2. * _geom->get_assembly_spacing())) * cm * .5);
0684   //
0685   //      G4LogicalVolume* wall_logic = new G4LogicalVolume(wall_solid, wall_mat,
0686   //                                                        G4String(G4String(GetName() + G4String("_SideWall_") + to_string(val.first))), 0, 0,
0687   //                                                        nullptr);
0688   //      GetDisplayAction()->AddVolume(wall_logic, "Wall");
0689   //
0690   //      const G4Transform3D wall_trans =
0691   //          G4TranslateZ3D(sign_z * (_geom->get_length() * cm / 4)) *
0692   //          G4TranslateY3D(enclosure_center) *
0693   //          G4TranslateX3D(sign_azimuth * enclosure_half_height_half_width) *
0694   //          G4RotateZ3D(azimuth_roate) *
0695   //          G4TranslateX3D(-sign_azimuth * (_geom->get_sidewall_thickness() * cm / 2.0 + _geom->get_sidewall_outer_torr() * cm));
0696   //
0697   //      G4PVPlacement* wall_phys = new G4PVPlacement(wall_trans, wall_logic,
0698   //                                                   G4String(GetName()) + G4String("_SideWall_") + to_string(val.first), sec_logic,
0699   //                                                   false, val.first, OverlapCheck());
0700   //
0701   //      calo_vol[wall_phys] = val.first;
0702   //      assert(gdml_config);
0703   //      gdml_config->exclude_physical_vol(wall_phys);
0704   //    }
0705   //  }
0706   //
0707   //  // construct dividers
0708   //  if (_geom->get_divider_width() > 0)
0709   //  {
0710   //    if (_geom->get_construction_verbose() >= 1)
0711   //    {
0712   //      cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg::" << GetName()
0713   //           << " - construct dividers" << endl;
0714   //    }
0715   //
0716   //    G4Material* divider_mat = G4Material::GetMaterial(_geom->get_divider_mat());
0717   //    assert(divider_mat);
0718   //
0719   //    int ID = 300;
0720   //    for (const auto& geom : divider_azimuth_geoms)
0721   //    {
0722   //      G4Box* divider_solid = new G4Box(G4String(GetName() + G4String("_Divider_") + to_string(ID)),
0723   //                                       geom.thickness / 2.0,
0724   //                                       geom.width / 2.,
0725   //                                       (_geom->get_length() / 2. - 2 * (_geom->get_sidewall_thickness() + 2. * _geom->get_assembly_spacing())) * cm * .5);
0726   //
0727   //      G4LogicalVolume* wall_logic = new G4LogicalVolume(divider_solid, divider_mat,
0728   //                                                        G4String(G4String(GetName() + G4String("_Divider_") + to_string(ID))), 0, 0,
0729   //                                                        nullptr);
0730   //      GetDisplayAction()->AddVolume(wall_logic, "Divider");
0731   //
0732   //      for (int sign_z = -1; sign_z <= 1; sign_z += 2)
0733   //      {
0734   //        G4Transform3D wall_trans =
0735   //            G4TranslateX3D(geom.projection_center_x) *
0736   //            G4TranslateY3D(geom.projection_center_y) *
0737   //            G4RotateZ3D(geom.angle) *
0738   //            G4TranslateY3D(geom.radial_displacement) *
0739   //            G4TranslateZ3D(sign_z * (_geom->get_length() * cm / 4));
0740   //
0741   //        G4PVPlacement* wall_phys = new G4PVPlacement(wall_trans, wall_logic,
0742   //                                                     G4String(GetName()) + G4String("_Divider_") + to_string(ID), sec_logic,
0743   //                                                     false, ID, OverlapCheck());
0744   //
0745   //        calo_vol[wall_phys] = ID;
0746   //        //        assert(gdml_config);
0747   //        //        gdml_config->exclude_physical_vol(wall_phys);
0748   //
0749   //        if (Verbosity())
0750   //        {
0751   //          cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg - placing divider " << wall_phys->GetName() << " copy ID " << ID << endl;
0752   //        }
0753   //
0754   //        ++ID;
0755   //      }
0756   //    }  //    for (const auto & geom : divider_azimuth_geoms)
0757   //  }
0758 
0759   //  // construct towers
0760   //
0761   BOOST_FOREACH (const SpacalGeom_t::tower_map_t::value_type& val, _geom->get_sector_tower_map())
0762   {
0763     SpacalGeom_t::geom_tower g_tower = val.second;
0764 
0765     const int tower_id = g_tower.id;
0766     const int tower_phi_id_in_sec = tower_id % 10;
0767     assert(tower_phi_id_in_sec >= 0);
0768     assert(tower_phi_id_in_sec < phi_bin_in_sec);
0769 
0770     const auto& block_azimuth_geom = block_azimuth_geoms.at(tower_phi_id_in_sec);
0771 
0772     G4LogicalVolume* LV_tower = Construct_Tower(g_tower);
0773 
0774     G4Transform3D block_trans =
0775         G4TranslateX3D(block_azimuth_geom.projection_center_x) *
0776         G4TranslateY3D(block_azimuth_geom.projection_center_y) *
0777         G4RotateZ3D(block_azimuth_geom.angle) *
0778         G4TranslateX3D(g_tower.centralX * cm) *
0779         G4TranslateY3D(g_tower.centralY * cm) *
0780         G4TranslateZ3D(g_tower.centralZ * cm) *
0781         G4RotateX3D(g_tower.pRotationAngleX * rad);
0782 
0783     const bool overlapcheck_block = OverlapCheck() and (_geom->get_construction_verbose() >= 2);
0784 
0785     G4PVPlacement* block_phys = new G4PVPlacement(block_trans, LV_tower,
0786                                                   G4String(GetName()) + G4String("_Tower_") + to_string(g_tower.id), sec_logic, false,
0787                                                   g_tower.id, overlapcheck_block);
0788     block_vol[block_phys] = g_tower.id;
0789     //    assert(gdml_config);
0790     //    gdml_config->exclude_physical_vol(block_phys);
0791 
0792     if (g_tower.LightguideHeight > 0)
0793     {
0794       // also build a light guide
0795 
0796       for (int ix = 0; ix < g_tower.NSubtowerX; ix++)
0797       //  int ix = 0;
0798       {
0799         for (int iy = 0; iy < g_tower.NSubtowerY; iy++)
0800         //        int iy = 0;
0801         {
0802           G4LogicalVolume* LV_lg = Construct_LightGuide(g_tower, ix,
0803                                                         iy);
0804 
0805           G4PVPlacement* lg_phys = new G4PVPlacement(block_trans, LV_lg, LV_lg->GetName(),
0806                                                      sec_logic, false, g_tower.id, overlapcheck_block);
0807 
0808           block_vol[lg_phys] = g_tower.id * 100 + ix * 10 + iy;
0809 
0810           //          assert(gdml_config);
0811           //          gdml_config->exclude_physical_vol(lg_phys);
0812         }
0813       }
0814     }
0815   }
0816 
0817   cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg::" << GetName()
0818        << " - constructed " << _geom->get_sector_tower_map().size()
0819        << " unique towers" << endl;
0820 
0821   return make_pair(sec_logic, G4Transform3D::Identity);
0822 }
0823 
0824 G4LogicalVolume*
0825 PHG4SpacalPrototype4Detector::Construct_Fiber(const G4double length,
0826                                               const string& id)
0827 {
0828   G4Tubs* fiber_solid = new G4Tubs(G4String(GetName() + string("_fiber") + id),
0829                                    0, _geom->get_fiber_outer_r() * cm, length / 2.0, 0, twopi);
0830 
0831   G4Material* clading_mat = G4Material::GetMaterial(
0832       _geom->get_fiber_clading_mat());
0833   assert(clading_mat);
0834 
0835   G4LogicalVolume* fiber_logic = new G4LogicalVolume(fiber_solid, clading_mat,
0836                                                      G4String(G4String(GetName() + string("_fiber") + id)), 0, 0,
0837                                                      clading_step_limits);
0838 
0839   {
0840     G4VisAttributes* VisAtt = new G4VisAttributes();
0841     PHG4Utils::SetColour(VisAtt, "G4_POLYSTYRENE");
0842     VisAtt->SetVisibility(_geom->is_virualize_fiber());
0843     VisAtt->SetForceSolid(_geom->is_virualize_fiber());
0844     fiber_logic->SetVisAttributes(VisAtt);
0845   }
0846 
0847   G4Tubs* core_solid = new G4Tubs(
0848       G4String(GetName() + string("_fiber_core") + id), 0,
0849       _geom->get_fiber_core_diameter() * cm / 2, length / 2.0, 0, twopi);
0850 
0851   G4Material* core_mat = G4Material::GetMaterial(_geom->get_fiber_core_mat());
0852   assert(core_mat);
0853 
0854   G4LogicalVolume* core_logic = new G4LogicalVolume(core_solid, core_mat,
0855                                                     G4String(G4String(GetName() + string("_fiber_core") + id)), 0, 0,
0856                                                     fiber_core_step_limits);
0857 
0858   {
0859     G4VisAttributes* VisAtt = new G4VisAttributes();
0860     PHG4Utils::SetColour(VisAtt, "G4_POLYSTYRENE");
0861     VisAtt->SetVisibility(false);
0862     VisAtt->SetForceSolid(false);
0863     core_logic->SetVisAttributes(VisAtt);
0864   }
0865 
0866   const bool overlapcheck_fiber = OverlapCheck() and (_geom->get_construction_verbose() >= 3);
0867   G4PVPlacement* core_physi = new G4PVPlacement(0, G4ThreeVector(), core_logic,
0868                                                 G4String(G4String(GetName() + string("_fiber_core") + id)), fiber_logic,
0869                                                 false, 0, overlapcheck_fiber);
0870   fiber_core_vol[core_physi] = 0;
0871 
0872   return fiber_logic;
0873 }
0874 
0875 void PHG4SpacalPrototype4Detector::Print(const std::string& /*what*/) const
0876 {
0877   assert(_geom);
0878 
0879   cout << "PHG4SpacalPrototype4Detector::Print::" << GetName()
0880        << " - Print Geometry:" << endl;
0881   _geom->Print();
0882 
0883   return;
0884 }
0885 
0886 //! Fully projective spacal with 2D tapered modules. To speed up construction, same-length fiber is used cross one tower
0887 int PHG4SpacalPrototype4Detector::Construct_Fibers_SameLengthFiberPerTower(
0888     const PHG4SpacalPrototype4Detector::SpacalGeom_t::geom_tower& g_tower,
0889     G4LogicalVolume* LV_tower)
0890 {
0891   assert(_geom);
0892 
0893   // construct fibers
0894 
0895   // first check out the fibers geometry
0896 
0897   typedef map<int, pair<G4Vector3D, G4Vector3D> > fiber_par_map;
0898   fiber_par_map fiber_par;
0899   G4double min_fiber_length = g_tower.pDz * cm * 4;
0900 
0901   G4Vector3D v_zshift = G4Vector3D(tan(g_tower.pTheta) * cos(g_tower.pPhi),
0902                                    tan(g_tower.pTheta) * sin(g_tower.pPhi), 1) *
0903                         g_tower.pDz;
0904   //  int fiber_ID = 0;
0905   for (int ix = 0; ix < g_tower.NFiberX; ix++)
0906   //  int ix = 0;
0907   {
0908     const double weighted_ix = static_cast<double>(ix) / (g_tower.NFiberX - 1.);
0909 
0910     const double weighted_pDx1 = (g_tower.pDx1 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0911     const double weighted_pDx2 = (g_tower.pDx2 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0912 
0913     const double weighted_pDx3 = (g_tower.pDx3 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0914     const double weighted_pDx4 = (g_tower.pDx4 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0915 
0916     for (int iy = 0; iy < g_tower.NFiberY; iy++)
0917     //        int iy = 0;
0918     {
0919       if ((ix + iy) % 2 == 1)
0920         continue;  // make a triangle pattern
0921 
0922       const double weighted_iy = static_cast<double>(iy) / (g_tower.NFiberY - 1.);
0923 
0924       const double weighted_pDy1 = (g_tower.pDy1 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
0925       const double weighted_pDy2 = (g_tower.pDy2 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
0926 
0927       const double weighted_pDx12 = weighted_pDx1 * (1 - weighted_iy) + weighted_pDx2 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp1);
0928       const double weighted_pDx34 = weighted_pDx3 * (1 - weighted_iy) + weighted_pDx4 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp2);
0929 
0930       G4Vector3D v1 = G4Vector3D(weighted_pDx12, weighted_pDy1, 0) - v_zshift;
0931       G4Vector3D v2 = G4Vector3D(weighted_pDx34, weighted_pDy2, 0) + v_zshift;
0932 
0933       G4Vector3D vector_fiber = (v2 - v1);
0934       vector_fiber *= (vector_fiber.mag() - _geom->get_fiber_outer_r()) / vector_fiber.mag();  // shrink by fiber boundary protection
0935       G4Vector3D center_fiber = (v2 + v1) / 2;
0936 
0937       // convert to Geant4 units
0938       vector_fiber *= cm;
0939       center_fiber *= cm;
0940 
0941       const int fiber_ID = g_tower.compose_fiber_id(ix, iy);
0942       fiber_par[fiber_ID] = make_pair(vector_fiber, center_fiber);
0943 
0944       const G4double fiber_length = vector_fiber.mag();
0945 
0946       min_fiber_length = min(fiber_length, min_fiber_length);
0947 
0948       //          ++fiber_ID;
0949     }
0950   }
0951 
0952   int fiber_count = 0;
0953 
0954   const G4double fiber_length = min_fiber_length;
0955   vector<G4double> fiber_cut;
0956 
0957   ostringstream ss;
0958   ss << string("_Tower") << g_tower.id;
0959   G4LogicalVolume* fiber_logic = Construct_Fiber(fiber_length, ss.str());
0960 
0961   BOOST_FOREACH (const fiber_par_map::value_type& val, fiber_par)
0962   {
0963     const int fiber_ID = val.first;
0964     G4Vector3D vector_fiber = val.second.first;
0965     G4Vector3D center_fiber = val.second.second;
0966     const G4double optimal_fiber_length = vector_fiber.mag();
0967 
0968     const G4Vector3D v1 = center_fiber - 0.5 * vector_fiber;
0969 
0970     // keep a statistics
0971     assert(optimal_fiber_length - fiber_length >= 0);
0972     fiber_cut.push_back(optimal_fiber_length - fiber_length);
0973 
0974     center_fiber += (fiber_length / optimal_fiber_length - 1) * 0.5 * vector_fiber;
0975     vector_fiber *= fiber_length / optimal_fiber_length;
0976 
0977     //      const G4Vector3D v1_new = center_fiber - 0.5 *vector_fiber;
0978 
0979     if (_geom->get_construction_verbose() >= 3)
0980       cout
0981           << "PHG4SpacalPrototype4Detector::Construct_Fibers_SameLengthFiberPerTower::"
0982           << GetName() << " - constructed fiber " << fiber_ID << ss.str()
0983           //
0984           << ", Length = " << optimal_fiber_length << "-"
0985           << (optimal_fiber_length - fiber_length) << "mm, "  //
0986           << "x = " << center_fiber.x() << "mm, "             //
0987           << "y = " << center_fiber.y() << "mm, "             //
0988           << "z = " << center_fiber.z() << "mm, "             //
0989           << "vx = " << vector_fiber.x() << "mm, "            //
0990           << "vy = " << vector_fiber.y() << "mm, "            //
0991           << "vz = " << vector_fiber.z() << "mm, "            //
0992           << endl;
0993 
0994     const G4double rotation_angle = G4Vector3D(0, 0, 1).angle(vector_fiber);
0995     const G4Vector3D rotation_axis =
0996         rotation_angle == 0 ? G4Vector3D(1, 0, 0) : G4Vector3D(0, 0, 1).cross(vector_fiber);
0997 
0998     G4Transform3D fiber_place(
0999         G4Translate3D(center_fiber.x(), center_fiber.y(), center_fiber.z()) * G4Rotate3D(rotation_angle, rotation_axis));
1000 
1001     ostringstream name;
1002     name << GetName() + string("_Tower") << g_tower.id << "_fiber"
1003          << ss.str();
1004 
1005     const bool overlapcheck_fiber = OverlapCheck() and (_geom->get_construction_verbose() >= 3);
1006     G4PVPlacement* fiber_physi = new G4PVPlacement(fiber_place, fiber_logic,
1007                                                    G4String(name.str()), LV_tower, false, fiber_ID,
1008                                                    overlapcheck_fiber);
1009     fiber_vol[fiber_physi] = fiber_ID;
1010 
1011     fiber_count++;
1012   }
1013 
1014   if (_geom->get_construction_verbose() >= 2)
1015     cout
1016         << "PHG4SpacalPrototype4Detector::Construct_Fibers_SameLengthFiberPerTower::"
1017         << GetName() << " - constructed tower ID " << g_tower.id << " with "
1018         << fiber_count << " fibers. Average fiber length cut = "
1019         << accumulate(fiber_cut.begin(), fiber_cut.end(), 0.0) / fiber_cut.size() << " mm" << endl;
1020 
1021   return fiber_count;
1022 }
1023 
1024 //! a block along z axis built with G4Trd that is slightly tapered in x dimension
1025 G4LogicalVolume*
1026 PHG4SpacalPrototype4Detector::Construct_Tower(
1027     const PHG4SpacalPrototype4Detector::SpacalGeom_t::geom_tower& g_tower)
1028 {
1029   assert(_geom);
1030 
1031   if (_geom->get_construction_verbose() >= 2)
1032   {
1033     cout << "PHG4SpacalPrototype4Detector::Construct_Tower::" << GetName()
1034          << " - constructed tower ID " << g_tower.id
1035          << " with geometry parameter: ";
1036     g_tower.identify(cout);
1037   }
1038 
1039   std::ostringstream sout;
1040   sout << "_" << g_tower.id;
1041   const G4String sTowerID(sout.str());
1042 
1043   //Processed PostionSeeds 1 from 1 1
1044 
1045   G4Trap* block_solid = new G4Trap(
1046       /*const G4String& pName*/ G4String(GetName()) + sTowerID,
1047       g_tower.pDz * cm,                                         // G4double pDz,
1048       g_tower.pTheta * rad, g_tower.pPhi * rad,                 // G4double pTheta, G4double pPhi,
1049       g_tower.pDy1 * cm, g_tower.pDx1 * cm, g_tower.pDx2 * cm,  // G4double pDy1, G4double pDx1, G4double pDx2,
1050       g_tower.pAlp1 * rad,                                      // G4double pAlp1,
1051       g_tower.pDy2 * cm, g_tower.pDx3 * cm, g_tower.pDx4 * cm,  // G4double pDy2, G4double pDx3, G4double pDx4,
1052       g_tower.pAlp2 * rad                                       // G4double pAlp2 //
1053   );
1054 
1055   G4Material* cylinder_mat = G4Material::GetMaterial(
1056       _geom->get_absorber_mat());
1057   assert(cylinder_mat);
1058 
1059   G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, cylinder_mat,
1060                                                      G4String(G4String(GetName()) + string("_Tower") + sTowerID), 0, 0,
1061                                                      step_limits);
1062 
1063   G4VisAttributes* VisAtt = new G4VisAttributes();
1064   //  PHG4Utils::SetColour(VisAtt, "W_Epoxy");
1065   VisAtt->SetColor(.3, .3, .3, .3);
1066   VisAtt->SetVisibility(
1067       _geom->is_azimuthal_seg_visible() or _geom->is_virualize_fiber());
1068   VisAtt->SetForceSolid(not _geom->is_virualize_fiber());
1069   block_logic->SetVisAttributes(VisAtt);
1070 
1071   // construct fibers
1072 
1073   int fiber_count = 0;
1074 
1075   fiber_count = Construct_Fibers_SameLengthFiberPerTower(g_tower, block_logic);
1076 
1077   if (_geom->get_construction_verbose() >= 2)
1078     cout << "PHG4SpacalPrototype4Detector::Construct_Tower::" << GetName()
1079          << " - constructed tower ID " << g_tower.id << " with " << fiber_count
1080          << " fibers using Construct_Fibers_SameLengthFiberPerTower" << endl;
1081 
1082   return block_logic;
1083 }
1084 
1085 G4LogicalVolume*
1086 PHG4SpacalPrototype4Detector::Construct_LightGuide(
1087     const PHG4SpacalPrototype4Detector::SpacalGeom_t::geom_tower& g_tower,
1088     const int index_x, const int index_y)
1089 {
1090   assert(_geom);
1091 
1092   std::ostringstream sout;
1093   sout << "_Lightguide_" << g_tower.id << "_" << index_x << "_" << index_y;
1094   const G4String sTowerID(sout.str());
1095 
1096   assert(g_tower.LightguideHeight > 0);
1097 
1098   // light guide parameters in PHENIX units
1099   const double weight_x1 = 1 - (double) index_y / g_tower.NSubtowerY;
1100   const double weight_x2 = 1 - (double) (index_y + 1) / g_tower.NSubtowerY;
1101   const double weight_xcenter = 1 - (double) (index_y + 0.5) / g_tower.NSubtowerY;
1102 
1103   assert(weight_x1 >= 0 and weight_x1 <= 1);
1104   assert(weight_x2 >= 0 and weight_x2 <= 1);
1105   assert(weight_xcenter >= 0 and weight_xcenter <= 1);
1106 
1107   const double lg_pDx1 = (g_tower.pDx1 * weight_x1  //
1108                           + g_tower.pDx2 * (1 - weight_x1)) /
1109                          g_tower.NSubtowerX;
1110   const double lg_pDx2 = (g_tower.pDx1 * weight_x2  //
1111                           + g_tower.pDx2 * (1 - weight_x2)) /
1112                          g_tower.NSubtowerX;
1113   const double lg_pDy1 = g_tower.pDy1 / g_tower.NSubtowerY;
1114   const double lg_Alp1 = atan(
1115       (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));
1116 
1117   const double shift_xcenter = (g_tower.pDx1 * weight_xcenter           //
1118                                 + g_tower.pDx2 * (1 - weight_xcenter))  //
1119                                *                                        //
1120                                (-g_tower.NSubtowerX + 1. + 2 * index_x) / (double) (g_tower.NSubtowerX);
1121   const double shift_ycenter = g_tower.pDy1  //
1122                                *             //
1123                                (-g_tower.NSubtowerY + 1. + 2 * index_y) / (double) (g_tower.NSubtowerY);
1124 
1125   G4VSolid* block_solid = new G4Trap(
1126       /*const G4String& pName*/ G4String(GetName()) + sTowerID,
1127       0.5 * g_tower.LightguideHeight * cm,  // G4double pDz,
1128       0 * rad, 0 * rad,                     // G4double pTheta, G4double pPhi,
1129       g_tower.LightguideTaperRatio * lg_pDy1 * cm,
1130       g_tower.LightguideTaperRatio * lg_pDx1 * cm,
1131       g_tower.LightguideTaperRatio * lg_pDx2 * cm,  // G4double pDy1, G4double pDx1, G4double pDx2,
1132       lg_Alp1 * rad,                                // G4double pAlp1,
1133       lg_pDy1 * cm, lg_pDx1 * cm, lg_pDx2 * cm,     // G4double pDy2, G4double pDx3, G4double pDx4,
1134       lg_Alp1 * rad                                 // G4double pAlp2 //
1135   );
1136 
1137   block_solid = new G4DisplacedSolid(G4String(GetName() + "_displaced"),
1138                                      block_solid, 0,                                           //
1139                                      G4ThreeVector(                                            //
1140                                          tan(g_tower.pTheta * rad) * cos(g_tower.pPhi * rad),  //
1141                                          tan(g_tower.pTheta * rad) * sin(g_tower.pPhi * rad),  //
1142                                          1) *                                                  // G4ThreeVector
1143                                              -(g_tower.pDz) *
1144                                              cm                                                       //
1145                                          + G4ThreeVector(shift_xcenter * cm, shift_ycenter * cm, 0)   // shit in subtower direction
1146                                          + G4ThreeVector(0, 0, -0.5 * g_tower.LightguideHeight * cm)  //shift in the light guide height
1147   );
1148 
1149   G4Material* cylinder_mat = G4Material::GetMaterial(
1150       g_tower.LightguideMaterial);
1151   assert(cylinder_mat);
1152 
1153   G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, cylinder_mat,
1154                                                      G4String(G4String(GetName()) + string("_Tower") + sTowerID), 0, 0,
1155                                                      step_limits);
1156 
1157   G4VisAttributes* VisAtt = new G4VisAttributes();
1158   PHG4Utils::SetColour(VisAtt, g_tower.LightguideMaterial);
1159   //  VisAtt->SetColor(.3, .3, .3, .3);
1160   VisAtt->SetVisibility(
1161       _geom->is_azimuthal_seg_visible() or _geom->is_virualize_fiber());
1162   VisAtt->SetForceSolid(not _geom->is_virualize_fiber());
1163   block_logic->SetVisAttributes(VisAtt);
1164 
1165   return block_logic;
1166 }
[ 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!