g4simulation/g4detectors/PHG4FullProjSpacalDetector.cc

0001
0002 /*!
0003  * \file ${file_name}
0004  * \brief
0005  * \author Jin Huang <jhuang@bnl.gov>
0006  * \version $$Revision: 1.3 $$
0007  * \date $$Date: 2015/02/10 15:39:26 $$
0008  */
0009 #include "PHG4FullProjSpacalDetector.h"
0010
0011 #include "PHG4SpacalDisplayAction.h"
0012
0013 #include <g4gdml/PHG4GDMLConfig.hh>
0014
0015 #include <phool/recoConsts.h>
0016
0017 #include <Geant4/G4Box.hh>
0018 #include <Geant4/G4Exception.hh>          // for G4Exception, G4ExceptionD
0019 #include <Geant4/G4ExceptionSeverity.hh>  // for FatalException
0020 #include <Geant4/G4LogicalVolume.hh>
0021 #include <Geant4/G4PVPlacement.hh>
0022 #include <Geant4/G4PhysicalConstants.hh>
0023 #include <Geant4/G4String.hh>  // for G4String
0024 #include <Geant4/G4SystemOfUnits.hh>
0025 #include <Geant4/G4Transform3D.hh>
0026 #include <Geant4/G4Trap.hh>
0027 #include <Geant4/G4Tubs.hh>
0028 #include <Geant4/G4Types.hh>  // for G4double
0029 #include <Geant4/G4Vector3D.hh>
0030
0031 #include <TSystem.h>
0032
0033 #include <algorithm>
0034 #include <cassert>
0035 #include <cmath>
0036 #include <iostream>  // for operator<<, basic_ostream
0037 #include <map>       // for allocator, map<>::value_type
0038 #include <numeric>   // std::accumulate
0039 #include <sstream>
0040 #include <string>  // std::string, std::to_string
0041 #include <vector>  // for vector
0042
0043 class G4Material;
0044 class PHCompositeNode;
0045
0046 //_______________________________________________________________
0047 // note this inactive thickness is ~1.5% of a radiation length
0048 PHG4FullProjSpacalDetector::PHG4FullProjSpacalDetector(PHG4Subsystem* subsys, PHCompositeNode* Node,
0049                                                        const std::string& dnam, PHParameters* parameters, const int lyr)
0050   : PHG4SpacalDetector(subsys, Node, dnam, parameters, lyr, false)
0051 {
0052   assert(_geom == nullptr);
0053
0054   _geom = new SpacalGeom_t();
0055   if (_geom == nullptr)
0056   {
0057     std::cout
0058         << "PHG4FullProjSpacalDetector::Constructor - Fatal Error - invalid geometry object!"
0059         << std::endl;
0060     gSystem->Exit(1);
0061   }
0062
0063   // this class loads Chris Cullen 2D spacal design July 2015 by default.
0064   //  this step is deprecated now
0065   //  get_geom_v3()->load_demo_sector_tower_map_2015_Chris_Cullen_2D_spacal();
0066
0067   assert(parameters);
0068   get_geom_v3()->ImportParameters(*parameters);
0069
0070   //  std::cout <<"PHG4FullProjSpacalDetector::Constructor -  get_geom_v3()->Print();"<<std::endl;
0071   //  get_geom_v3()->Print();
0072 }
0073
0074 //_______________________________________________________________
0075 void PHG4FullProjSpacalDetector::ConstructMe(G4LogicalVolume* logicWorld)
0076 {
0077   if (get_geom_v3()->get_construction_verbose() >= 1)
0078   {
0079     std::cout << "PHG4FullProjSpacalDetector::Construct::" << GetName()
0080               << " - start with PHG4SpacalDetector::Construct()." << std::endl;
0081   }
0082
0083   PHG4SpacalDetector::ConstructMe(logicWorld);
0084
0085   if (get_geom_v3()->get_construction_verbose() >= 1)
0086   {
0087     std::cout << "PHG4FullProjSpacalDetector::Construct::" << GetName()
0088               << " - Completed." << std::endl;
0089   }
0090 }
0091
0092 std::pair<G4LogicalVolume*, G4Transform3D>
0093 PHG4FullProjSpacalDetector::Construct_AzimuthalSeg()
0094 {
0095   if (!(get_geom_v3()->get_azimuthal_n_sec() > 4))
0096   {
0097     std::cout << "azimuthal_n_sec <= 4: " << get_geom_v3()->get_azimuthal_n_sec() << std::endl;
0098     gSystem->Exit(1);
0099   }
0100
0101   G4Tubs* sec_solid = new G4Tubs(G4String(GetName() + std::string("_sec")),
0102                                  get_geom_v3()->get_radius() * cm, get_geom_v3()->get_max_radius() * cm,
0103                                  get_geom_v3()->get_length() * cm / 2.0,
0104                                  halfpi - pi / get_geom_v3()->get_azimuthal_n_sec(),
0105                                  twopi / get_geom_v3()->get_azimuthal_n_sec());
0106
0107   recoConsts* rc = recoConsts::instance();
0108   G4Material* cylinder_mat = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
0109   assert(cylinder_mat);
0110
0111   G4LogicalVolume* sec_logic = new G4LogicalVolume(sec_solid, cylinder_mat,
0112                                                    G4String(G4String(GetName() + std::string("_sec"))), nullptr, nullptr);
0113
0114   GetDisplayAction()->AddVolume(sec_logic, "Sector");
0115
0116   // construct walls
0117
0118   G4Material* wall_mat = GetDetectorMaterial(get_geom_v3()->get_sidewall_mat());
0119   assert(wall_mat);
0120
0121   if (get_geom_v3()->get_sidewall_thickness() > 0)
0122   {
0123     // end walls
0124     if (get_geom_v3()->get_construction_verbose() >= 1)
0125     {
0126       std::cout << "PHG4FullProjSpacalDetector::Construct_AzimuthalSeg::" << GetName()
0127                 << " - construct end walls." << std::endl;
0128     }
0129     G4Tubs* wall_solid = new G4Tubs(G4String(GetName() + std::string("_EndWall")),
0130                                     get_geom_v3()->get_radius() * cm + get_geom_v3()->get_sidewall_outer_torr() * cm,
0131                                     get_geom_v3()->get_max_radius() * cm - get_geom_v3()->get_sidewall_outer_torr() * cm,
0132                                     get_geom_v3()->get_sidewall_thickness() * cm / 2.0,
0133                                     halfpi - pi / get_geom_v3()->get_azimuthal_n_sec(),
0134                                     twopi / get_geom_v3()->get_azimuthal_n_sec());
0135
0136     G4LogicalVolume* wall_logic = new G4LogicalVolume(wall_solid, wall_mat,
0137                                                       G4String(G4String(GetName() + std::string("_EndWall"))), nullptr, nullptr,
0138                                                       nullptr);
0139     GetDisplayAction()->AddVolume(wall_logic, "WallProj");
0140
0141     using z_locations_t = std::map<int, double>;
0142     z_locations_t z_locations;
0143     z_locations[1000] = get_geom_v3()->get_sidewall_thickness() * cm / 2.0 + get_geom_v3()->get_assembly_spacing() * cm;
0144     z_locations[1001] = get_geom_v3()->get_length() * cm / 2.0 - (get_geom_v3()->get_sidewall_thickness() * cm / 2.0 + get_geom_v3()->get_assembly_spacing() * cm);
0145     z_locations[1100] = -(get_geom_v3()->get_sidewall_thickness() * cm / 2.0 + get_geom_v3()->get_assembly_spacing() * cm);
0146     z_locations[1101] = -(get_geom_v3()->get_length() * cm / 2.0 - (get_geom_v3()->get_sidewall_thickness() * cm / 2.0 + get_geom_v3()->get_assembly_spacing() * cm));
0147
0148     for (z_locations_t::value_type& val : z_locations)
0149     {
0150       if (get_geom_v3()->get_construction_verbose() >= 2)
0151       {
0152         std::cout << "PHG4FullProjSpacalDetector::Construct_AzimuthalSeg::"
0153                   << GetName() << " - constructed End Wall ID " << val.first
0154                   << " @ Z = " << val.second << std::endl;
0155       }
0156
0157       G4Transform3D wall_trans = G4TranslateZ3D(val.second);
0158
0159       G4PVPlacement* wall_phys = new G4PVPlacement(wall_trans, wall_logic,
0160                                                    G4String(GetName()) + G4String("_EndWall"), sec_logic,
0161                                                    false, val.first, OverlapCheck());
0162
0163       calo_vol[wall_phys] = val.first;
0164       assert(gdml_config);
0165       gdml_config->exclude_physical_vol(wall_phys);
0166     }
0167   }
0168
0169   if (get_geom_v3()->get_sidewall_thickness() > 0)
0170   {
0171     // side walls
0172     if (get_geom_v3()->get_construction_verbose() >= 1)
0173     {
0174       std::cout << "PHG4FullProjSpacalDetector::Construct_AzimuthalSeg::" << GetName()
0175                 << " - construct side walls." << std::endl;
0176     }
0177     G4Box* wall_solid = new G4Box(G4String(GetName() + std::string("_SideWall")),
0178                                   get_geom_v3()->get_sidewall_thickness() * cm / 2.0,
0179                                   get_geom_v3()->get_thickness() * cm / 2. - 2 * get_geom_v3()->get_sidewall_outer_torr() * cm,
0180                                   (get_geom_v3()->get_length() / 2. - 2 * (get_geom_v3()->get_sidewall_thickness() + 2. * get_geom_v3()->get_assembly_spacing())) * cm * .5);
0181
0182     G4LogicalVolume* wall_logic = new G4LogicalVolume(wall_solid, wall_mat,
0183                                                       G4String(G4String(GetName() + std::string("_SideWall"))), nullptr, nullptr,
0184                                                       nullptr);
0185     GetDisplayAction()->AddVolume(wall_logic, "WallProj");
0186
0187     using sign_t = std::map<int, std::pair<int, int>>;
0188     sign_t signs;
0189     signs[2000] = std::make_pair(+1, +1);
0190     signs[2001] = std::make_pair(+1, -1);
0191     signs[2100] = std::make_pair(-1, +1);
0192     signs[2101] = std::make_pair(-1, -1);
0193
0194     for (sign_t::value_type& val : signs)
0195     {
0196       const int sign_z = val.second.first;
0197       const int sign_azimuth = val.second.second;
0198
0199       if (get_geom_v3()->get_construction_verbose() >= 2)
0200       {
0201         std::cout << "PHG4FullProjSpacalDetector::Construct_AzimuthalSeg::"
0202                   << GetName() << " - constructed Side Wall ID " << val.first
0203                   << " with"
0204                   << " Shift X = "
0205                   << sign_azimuth * (get_geom_v3()->get_sidewall_thickness() * cm / 2.0 + get_geom_v3()->get_sidewall_outer_torr() * cm)
0206                   << " Rotation Z = "
0207                   << sign_azimuth * pi / get_geom_v3()->get_azimuthal_n_sec()
0208                   << " Shift Z = " << sign_z * (get_geom_v3()->get_length() * cm / 4)
0209                   << std::endl;
0210       }
0211
0212       G4Transform3D wall_trans = G4RotateZ3D(
0213                                      sign_azimuth * pi / get_geom_v3()->get_azimuthal_n_sec()) *
0214                                  G4TranslateZ3D(sign_z * (get_geom_v3()->get_length() * cm / 4)) * G4TranslateY3D(get_geom_v3()->get_radius() * cm + get_geom_v3()->get_thickness() * cm / 2.) * G4TranslateX3D(sign_azimuth * (get_geom_v3()->get_sidewall_thickness() * cm / 2.0 + get_geom_v3()->get_sidewall_outer_torr() * cm));
0215
0216       G4PVPlacement* wall_phys = new G4PVPlacement(wall_trans, wall_logic,
0217                                                    G4String(GetName()) + G4String("_EndWall"), sec_logic,
0218                                                    false, val.first, OverlapCheck());
0219
0220       calo_vol[wall_phys] = val.first;
0221
0222       assert(gdml_config);
0223       gdml_config->exclude_physical_vol(wall_phys);
0224     }
0225   }
0226
0227   // construct towers
0228
0229   for (const SpacalGeom_t::tower_map_t::value_type& val : get_geom_v3()->get_sector_tower_map())
0230   {
0231     const SpacalGeom_t::geom_tower& g_tower = val.second;
0232     G4LogicalVolume* LV_tower = Construct_Tower(g_tower);
0233
0234     G4Transform3D block_trans = G4TranslateX3D(g_tower.centralX * cm) * G4TranslateY3D(g_tower.centralY * cm) * G4TranslateZ3D(g_tower.centralZ * cm) * G4RotateX3D(g_tower.pRotationAngleX * rad);
0235
0236     const bool overlapcheck_block = OverlapCheck() and (get_geom_v3()->get_construction_verbose() >= 2);
0237
0238     G4PVPlacement* block_phys = new G4PVPlacement(block_trans, LV_tower,
0239                                                   G4String(GetName()) + G4String("_Tower"), sec_logic, false,
0240                                                   g_tower.id, overlapcheck_block);
0241     block_vol[block_phys] = g_tower.id;
0242
0243     assert(gdml_config);
0244     gdml_config->exclude_physical_vol(block_phys);
0245   }
0246
0247   std::cout << "PHG4FullProjSpacalDetector::Construct_AzimuthalSeg::" << GetName()
0248             << " - constructed " << get_geom_v3()->get_sector_tower_map().size()
0249             << " unique towers" << std::endl;
0250
0251   return std::make_pair(sec_logic, G4Transform3D::Identity);
0252 }
0253
0254 //! Fully projective spacal with 2D tapered modules. To speed up construction, same-length fiber is used cross one tower
0255 int PHG4FullProjSpacalDetector::Construct_Fibers_SameLengthFiberPerTower(
0256     const PHG4FullProjSpacalDetector::SpacalGeom_t::geom_tower& g_tower,
0257     G4LogicalVolume* LV_tower)
0258 {
0259   // construct fibers
0260
0261   // first check out the fibers geometry
0262
0263   using fiber_par_map = std::map<int, std::pair<G4Vector3D, G4Vector3D>>;
0264   fiber_par_map fiber_par;
0265   G4double min_fiber_length = g_tower.pDz * cm * 4;
0266
0267   G4Vector3D v_zshift = G4Vector3D(tan(g_tower.pTheta) * cos(g_tower.pPhi),
0268                                    tan(g_tower.pTheta) * sin(g_tower.pPhi), 1) *
0269                         g_tower.pDz;
0270   //  int fiber_ID = 0;
0271   for (int ix = 0; ix < g_tower.NFiberX; ix++)
0272   //  int ix = 0;
0273   {
0274     const double weighted_ix = static_cast<double>(ix) / (g_tower.NFiberX - 1.);
0275
0276     const double weighted_pDx1 = (g_tower.pDx1 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0277     const double weighted_pDx2 = (g_tower.pDx2 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0278
0279     const double weighted_pDx3 = (g_tower.pDx3 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0280     const double weighted_pDx4 = (g_tower.pDx4 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0281
0282     for (int iy = 0; iy < g_tower.NFiberY; iy++)
0283     //        int iy = 0;
0284     {
0285       if ((ix + iy) % 2 == 1)
0286       {
0287         continue;  // make a triangle pattern
0288       }
0289
0290       const double weighted_iy = static_cast<double>(iy) / (g_tower.NFiberY - 1.);
0291
0292       const double weighted_pDy1 = (g_tower.pDy1 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
0293       const double weighted_pDy2 = (g_tower.pDy2 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
0294
0295       const double weighted_pDx12 = weighted_pDx1 * (1 - weighted_iy) + weighted_pDx2 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp1);
0296       const double weighted_pDx34 = weighted_pDx3 * (1 - weighted_iy) + weighted_pDx4 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp2);
0297
0298       G4Vector3D v1 = G4Vector3D(weighted_pDx12, weighted_pDy1, 0) - v_zshift;
0299       G4Vector3D v2 = G4Vector3D(weighted_pDx34, weighted_pDy2, 0) + v_zshift;
0300
0301       G4Vector3D vector_fiber = (v2 - v1);
0302       vector_fiber *= (vector_fiber.mag() - get_geom_v3()->get_fiber_outer_r()) / vector_fiber.mag();  // shrink by fiber boundary protection
0303       G4Vector3D center_fiber = (v2 + v1) / 2;
0304
0305       // convert to Geant4 units
0306       vector_fiber *= cm;
0307       center_fiber *= cm;
0308
0309       const int fiber_ID = g_tower.compose_fiber_id(ix, iy);
0310       fiber_par[fiber_ID] = std::make_pair(vector_fiber,
0311                                            center_fiber);
0312
0313       const G4double fiber_length = vector_fiber.mag();
0314
0315       min_fiber_length = std::min(fiber_length, min_fiber_length);
0316
0317       //          ++fiber_ID;
0318     }
0319   }
0320
0321   int fiber_count = 0;
0322
0323   const G4double fiber_length = min_fiber_length;
0324   std::vector<G4double> fiber_cut;
0325
0326   std::stringstream ss;
0327   ss << std::string("_Tower") << g_tower.id;
0328   G4LogicalVolume* fiber_logic = Construct_Fiber(fiber_length, ss.str());
0329
0330   for (const fiber_par_map::value_type& val : fiber_par)
0331   {
0332     const int fiber_ID = val.first;
0333     G4Vector3D vector_fiber = val.second.first;
0334     G4Vector3D center_fiber = val.second.second;
0335     const G4double optimal_fiber_length = vector_fiber.mag();
0336
0337     const G4Vector3D v1 = center_fiber - 0.5 * vector_fiber;
0338
0339     // keep a statistics
0340     assert(optimal_fiber_length - fiber_length >= 0);
0341     fiber_cut.push_back(optimal_fiber_length - fiber_length);
0342
0343     center_fiber += (fiber_length / optimal_fiber_length - 1) * 0.5 * vector_fiber;
0344     vector_fiber *= fiber_length / optimal_fiber_length;
0345
0346     //      const G4Vector3D v1_new = center_fiber - 0.5 *vector_fiber;
0347
0348     if (get_geom_v3()->get_construction_verbose() >= 3)
0349     {
0350       std::cout << "PHG4FullProjSpacalDetector::Construct_Fibers_SameLengthFiberPerTower::" << GetName()
0351                 << " - constructed fiber " << fiber_ID << ss.str()  //
0352                 << ", Length = " << optimal_fiber_length << "-"
0353                 << (optimal_fiber_length - fiber_length) << "mm, "  //
0354                 << "x = " << center_fiber.x() << "mm, "             //
0355                 << "y = " << center_fiber.y() << "mm, "             //
0356                 << "z = " << center_fiber.z() << "mm, "             //
0357                 << "vx = " << vector_fiber.x() << "mm, "            //
0358                 << "vy = " << vector_fiber.y() << "mm, "            //
0359                 << "vz = " << vector_fiber.z() << "mm, "            //
0360                 << std::endl;
0361     }
0362
0363     const G4double rotation_angle = G4Vector3D(0, 0, 1).angle(vector_fiber);
0364     const G4Vector3D rotation_axis =
0365         rotation_angle == 0 ? G4Vector3D(1, 0, 0) : G4Vector3D(0, 0, 1).cross(vector_fiber);
0366
0367     G4Transform3D fiber_place(
0368         G4Translate3D(center_fiber.x(), center_fiber.y(), center_fiber.z()) * G4Rotate3D(rotation_angle, rotation_axis));
0369
0370     std::stringstream name;
0371     name << GetName() + std::string("_Tower") << g_tower.id << "_fiber"
0372          << ss.str();
0373
0374     const bool overlapcheck_fiber = OverlapCheck() and (get_geom_v3()->get_construction_verbose() >= 3);
0375     G4PVPlacement* fiber_physi = new G4PVPlacement(fiber_place, fiber_logic,
0376                                                    G4String(name.str()), LV_tower, false, fiber_ID,
0377                                                    overlapcheck_fiber);
0378     fiber_vol[fiber_physi] = fiber_ID;
0379
0380     assert(gdml_config);
0381     gdml_config->exclude_physical_vol(fiber_physi);
0382
0383     fiber_count++;
0384   }
0385
0386   if (get_geom_v3()->get_construction_verbose() >= 2)
0387   {
0388     std::cout
0389         << "PHG4FullProjSpacalDetector::Construct_Fibers_SameLengthFiberPerTower::"
0390         << GetName() << " - constructed tower ID " << g_tower.id << " with "
0391         << fiber_count << " fibers. Average fiber length cut = "
0392         << accumulate(fiber_cut.begin(), fiber_cut.end(), 0.0) / fiber_cut.size() << " mm" << std::endl;
0393   }
0394
0395   return fiber_count;
0396 }
0397
0398 //! a block along z axis built with G4Trd that is slightly tapered in x dimension
0399 int PHG4FullProjSpacalDetector::Construct_Fibers(
0400     const PHG4FullProjSpacalDetector::SpacalGeom_t::geom_tower& g_tower,
0401     G4LogicalVolume* LV_tower)
0402 {
0403   G4Vector3D v_zshift = G4Vector3D(tan(g_tower.pTheta) * cos(g_tower.pPhi),
0404                                    tan(g_tower.pTheta) * sin(g_tower.pPhi), 1) *
0405                         g_tower.pDz;
0406   int fiber_cnt = 0;
0407   for (int ix = 0; ix < g_tower.NFiberX; ix++)
0408   {
0409     const double weighted_ix = static_cast<double>(ix) / (g_tower.NFiberX - 1.);
0410
0411     const double weighted_pDx1 = (g_tower.pDx1 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0412     const double weighted_pDx2 = (g_tower.pDx2 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0413
0414     const double weighted_pDx3 = (g_tower.pDx3 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0415     const double weighted_pDx4 = (g_tower.pDx4 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0416
0417     for (int iy = 0; iy < g_tower.NFiberY; iy++)
0418     {
0419       if ((ix + iy) % 2 == 1)
0420       {
0421         continue;  // make a triangle pattern
0422       }
0423       const int fiber_ID = g_tower.compose_fiber_id(ix, iy);
0424
0425       const double weighted_iy = static_cast<double>(iy) / (g_tower.NFiberY - 1.);
0426
0427       const double weighted_pDy1 = (g_tower.pDy1 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
0428       const double weighted_pDy2 = (g_tower.pDy2 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
0429
0430       const double weighted_pDx12 = weighted_pDx1 * (1 - weighted_iy) + weighted_pDx2 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp1);
0431       const double weighted_pDx34 = weighted_pDx3 * (1 - weighted_iy) + weighted_pDx4 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp2);
0432
0433       G4Vector3D v1 = G4Vector3D(weighted_pDx12, weighted_pDy1, 0) - v_zshift;
0434       G4Vector3D v2 = G4Vector3D(weighted_pDx34, weighted_pDy2, 0) + v_zshift;
0435
0436       G4Vector3D vector_fiber = (v2 - v1);
0437       vector_fiber *= (vector_fiber.mag() - get_geom_v3()->get_fiber_outer_r()) / vector_fiber.mag();  // shrink by fiber boundary protection
0438       G4Vector3D center_fiber = (v2 + v1) / 2;
0439
0440       // convert to Geant4 units
0441       vector_fiber *= cm;
0442       center_fiber *= cm;
0443
0444       const G4double fiber_length = vector_fiber.mag();
0445
0446       std::stringstream ss;
0447       ss << std::string("_Tower") << g_tower.id;
0448       ss << "_x" << ix;
0449       ss << "_y" << iy;
0450       G4LogicalVolume* fiber_logic = Construct_Fiber(fiber_length,
0451                                                      ss.str());
0452
0453       if (get_geom_v3()->get_construction_verbose() >= 3)
0454       {
0455         std::cout << "PHG4FullProjSpacalDetector::Construct_Fibers::" << GetName()
0456                   << " - constructed fiber " << fiber_ID << ss.str()  //
0457                   << ", Length = " << fiber_length << "mm, "          //
0458                   << "x = " << center_fiber.x() << "mm, "             //
0459                   << "y = " << center_fiber.y() << "mm, "             //
0460                   << "z = " << center_fiber.z() << "mm, "             //
0461                   << "vx = " << vector_fiber.x() << "mm, "            //
0462                   << "vy = " << vector_fiber.y() << "mm, "            //
0463                   << "vz = " << vector_fiber.z() << "mm, "            //
0464                   << std::endl;
0465       }
0466
0467       const G4double rotation_angle = G4Vector3D(0, 0, 1).angle(
0468           vector_fiber);
0469       const G4Vector3D rotation_axis =
0470           rotation_angle == 0 ? G4Vector3D(1, 0, 0) : G4Vector3D(0, 0, 1).cross(vector_fiber);
0471
0472       G4Transform3D fiber_place(
0473           G4Translate3D(center_fiber.x(), center_fiber.y(),
0474                         center_fiber.z()) *
0475           G4Rotate3D(rotation_angle, rotation_axis));
0476
0477       std::stringstream name;
0478       name << GetName() + std::string("_Tower") << g_tower.id << "_fiber"
0479            << ss.str();
0480
0481       const bool overlapcheck_fiber = OverlapCheck() and (get_geom_v3()->get_construction_verbose() >= 3);
0482       G4PVPlacement* fiber_physi = new G4PVPlacement(fiber_place,
0483                                                      fiber_logic, G4String(name.str()), LV_tower, false,
0484                                                      fiber_ID, overlapcheck_fiber);
0485       fiber_vol[fiber_physi] = fiber_ID;
0486
0487       assert(gdml_config);
0488       gdml_config->exclude_physical_vol(fiber_physi);
0489
0490       ++fiber_cnt;
0491     }
0492   }
0493
0494   if (get_geom_v3()->get_construction_verbose() >= 3)
0495   {
0496     std::cout << "PHG4FullProjSpacalDetector::Construct_Fibers::" << GetName()
0497               << " - constructed tower ID " << g_tower.id << " with " << fiber_cnt
0498               << " fibers" << std::endl;
0499   }
0500
0501   return fiber_cnt;
0502 }
0503
0504 //! a block along z axis built with G4Trd that is slightly tapered in x dimension
0505 G4LogicalVolume*
0506 PHG4FullProjSpacalDetector::Construct_Tower(
0507     const PHG4FullProjSpacalDetector::SpacalGeom_t::geom_tower& g_tower)
0508 {
0509   std::stringstream sout;
0510   sout << "_" << g_tower.id;
0511   const G4String sTowerID(sout.str());
0512
0513   // Processed PostionSeeds 1 from 1 1
0514
0515   G4Trap* block_solid = new G4Trap(
0516       /*const G4String& pName*/ G4String(GetName()) + sTowerID,
0517       g_tower.pDz * cm,                                         // G4double pDz,
0518       g_tower.pTheta * rad, g_tower.pPhi * rad,                 // G4double pTheta, G4double pPhi,
0519       g_tower.pDy1 * cm, g_tower.pDx1 * cm, g_tower.pDx2 * cm,  // G4double pDy1, G4double pDx1, G4double pDx2,
0520       g_tower.pAlp1 * rad,                                      // G4double pAlp1,
0521       g_tower.pDy2 * cm, g_tower.pDx3 * cm, g_tower.pDx4 * cm,  // G4double pDy2, G4double pDx3, G4double pDx4,
0522       g_tower.pAlp2 * rad                                       // G4double pAlp2 //
0523   );
0524
0525   G4Material* cylinder_mat = GetDetectorMaterial(get_geom_v3()->get_absorber_mat());
0526   assert(cylinder_mat);
0527
0528   G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, cylinder_mat,
0529                                                      G4String(G4String(GetName()) + std::string("_Tower") + sTowerID), nullptr, nullptr,
0530                                                      nullptr);
0531
0532   GetDisplayAction()->AddVolume(block_logic, "Block");
0533
0534   // construct fibers
0535
0536   if (get_geom_v3()->get_config() == SpacalGeom_t::kFullProjective_2DTaper)
0537   {
0538     int fiber_count = Construct_Fibers(g_tower, block_logic);
0539
0540     if (get_geom_v3()->get_construction_verbose() >= 2)
0541     {
0542       std::cout << "PHG4FullProjSpacalDetector::Construct_Tower::" << GetName()
0543                 << " - constructed tower ID " << g_tower.id << " with "
0544                 << fiber_count << " fibers using Construct_Fibers" << std::endl;
0545     }
0546   }
0547   else if (get_geom_v3()->get_config() == SpacalGeom_t::kFullProjective_2DTaper_SameLengthFiberPerTower)
0548   {
0549     int fiber_count = Construct_Fibers_SameLengthFiberPerTower(g_tower,
0550                                                                block_logic);
0551
0552     if (get_geom_v3()->get_construction_verbose() >= 2)
0553     {
0554       std::cout << "PHG4FullProjSpacalDetector::Construct_Tower::" << GetName()
0555                 << " - constructed tower ID " << g_tower.id << " with "
0556                 << fiber_count
0557                 << " fibers using Construct_Fibers_SameLengthFiberPerTower" << std::endl;
0558     }
0559   }
0560   else
0561   {
0562     G4ExceptionDescription message;
0563     message << "can not recognize configuration type " << get_geom_v3()->get_config();
0564
0565     G4Exception("PHG4FullProjSpacalDetector::Construct_Tower", "Wrong",
0566                 FatalException, message, "");
0567   }
0568
0569   return block_logic;
0570 }
0571
0572 void PHG4FullProjSpacalDetector::Print(const std::string& /*what*/) const
0573 {
0574   std::cout << "PHG4FullProjSpacalDetector::Print::" << GetName()
0575             << " - Print Geometry:" << std::endl;
0576   get_geom_v3()->Print();
0577
0578   return;
0579 }