Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4FullProjTiltedSpacalDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4FullProjTiltedSpacalDetector.cc
1 // $$Id: PHG4FullProjTiltedSpacalDetector.cc,v 1.3 2015/02/10 15:39:26 pinkenbu Exp $$
2 
10 #include "PHG4FullProjTiltedSpacalDetector.h"
11 
12 #include "PHG4SpacalDisplayAction.h"
13 
14 #include <g4gdml/PHG4GDMLConfig.hh>
15 
16 #include <phool/recoConsts.h>
17 
18 #include <phparameter/PHParameters.h>
19 
20 #include <Geant4/G4Box.hh>
21 #include <Geant4/G4DisplacedSolid.hh>
22 #include <Geant4/G4Exception.hh> // for G4Exception
23 #include <Geant4/G4ExceptionSeverity.hh> // for FatalException
24 #include <Geant4/G4LogicalVolume.hh>
25 #include <Geant4/G4PVPlacement.hh>
26 #include <Geant4/G4PhysicalConstants.hh>
27 #include <Geant4/G4String.hh> // for G4String
28 #include <Geant4/G4SystemOfUnits.hh>
29 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
30 #include <Geant4/G4Transform3D.hh> // for G4Transform3D, G4TranslateY3D
31 #include <Geant4/G4Trap.hh>
32 #include <Geant4/G4Types.hh> // for G4double
33 #include <Geant4/G4Vector3D.hh>
34 
35 #include <TSystem.h>
36 
37 #include <algorithm>
38 #include <cassert>
39 #include <cmath>
40 #include <iostream> // for operator<<, basic_ostream
41 #include <limits> // for numeric_limits
42 #include <map> // for map<>::value_type, map
43 #include <memory> // for allocator_traits<>::value_...
44 #include <numeric> // std::accumulate
45 #include <sstream>
46 #include <string> // std::string, std::to_string
47 #include <vector> // for vector
48 
49 class G4Material;
50 class G4VSolid;
51 class PHCompositeNode;
52 
53 //_______________________________________________________________
54 // note this inactive thickness is ~1.5% of a radiation length
55 PHG4FullProjTiltedSpacalDetector::PHG4FullProjTiltedSpacalDetector(PHG4Subsystem* subsys, PHCompositeNode* Node,
56  const std::string& dnam, PHParameters* parameters, const int lyr)
57  : PHG4SpacalDetector(subsys, Node, dnam, parameters, lyr, false)
58  , m_Params(parameters)
59 {
60  assert(_geom == nullptr);
61 
62  _geom = new SpacalGeom_t();
63  if (_geom == nullptr)
64  {
65  std::cout
66  << "PHG4FullProjTiltedSpacalDetector::Constructor - Fatal Error - invalid geometry object!"
67  << std::endl;
68  gSystem->Exit(1);
69  }
70  assert(parameters);
71 
72  assert(parameters);
73  get_geom_v3()->ImportParameters(*parameters);
74  // std::cout <<"PHG4FullProjTiltedSpacalDetector::Constructor - get_geom_v3()->Print();"<<std::endl;
75  // get_geom_v3()->Print();
76 }
77 
78 //_______________________________________________________________
79 void PHG4FullProjTiltedSpacalDetector::ConstructMe(G4LogicalVolume* logicWorld)
80 {
81  if (get_geom_v3()->get_construction_verbose() >= 1)
82  {
83  std::cout << "PHG4FullProjTiltedSpacalDetector::Construct::" << GetName()
84  << " - start with PHG4SpacalDetector::Construct()." << std::endl;
85  }
86 
87  PHG4SpacalDetector::ConstructMe(logicWorld);
88 
89  if (get_geom_v3()->get_construction_verbose() >= 1)
90  {
91  std::cout << "PHG4FullProjTiltedSpacalDetector::Construct::" << GetName()
92  << " - Completed." << std::endl;
93  }
94 
95  if (m_Params->get_int_param("saveg4hit")) return;
96  try
97  {
98  AddCellGeometryNode();
99  AddTowerGeometryNode();
100  }
101  catch (std::exception& e)
102  {
103  std::cout << e.what() << std::endl;
104  // exit(1);
105  }
106 }
107 
108 std::pair<G4LogicalVolume*, G4Transform3D>
109 PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg()
110 {
111  if (!(get_geom_v3()->get_azimuthal_n_sec() > 4))
112  {
113  std::cout << "azimuthal n sec <= 4: " << get_geom_v3()->get_azimuthal_n_sec() << std::endl;
114  gSystem->Exit(1);
115  }
116 
117  // basic tilt geometry
118  const G4double half_chord_backend =
119  get_geom_v3()->get_max_radius() * cm * tan(pi / get_geom_v3()->get_azimuthal_n_sec()) //
120  + fabs(get_geom_v3()->get_thickness() * cm * 0.5 * tan(get_geom_v3()->get_azimuthal_tilt()));
121  const G4double reduced_outer_radius = sqrt(pow(get_geom_v3()->get_max_radius() * cm, 2) - half_chord_backend * half_chord_backend);
122  const G4double enclosure_depth = reduced_outer_radius - get_geom_v3()->get_radius() * cm;
123  const G4double enclosure_center = 0.5 * (reduced_outer_radius + get_geom_v3()->get_radius() * cm);
124  const G4double enclosure_half_height_half_width = enclosure_center * tan(pi / get_geom_v3()->get_azimuthal_n_sec());
125 
126  const G4double width_adj1 = tan(get_geom_v3()->get_azimuthal_tilt() - pi / get_geom_v3()->get_azimuthal_n_sec()) * enclosure_depth * 0.5;
127  const G4double width_adj2 = tan(get_geom_v3()->get_azimuthal_tilt() + pi / get_geom_v3()->get_azimuthal_n_sec()) * enclosure_depth * 0.5;
128 
129  const G4double center_adj = (width_adj1 + width_adj2) * 0.5;
130  const G4double center_tilt_angle = atan2(center_adj, enclosure_depth * 0.5);
131  const G4double inner_half_width = enclosure_half_height_half_width + 0.5 * (width_adj1 - width_adj2);
132  const G4double outter_half_width = enclosure_half_height_half_width + 0.5 * (-width_adj1 + width_adj2);
133 
134  // enclosure walls
135  const G4double edge1_tilt_angle = atan2(width_adj1, enclosure_depth * 0.5);
136  const G4double edge2_tilt_angle = atan2(width_adj2, enclosure_depth * 0.5);
137  const G4double edge1_half_depth = sqrt(width_adj1 * width_adj1 + enclosure_depth * enclosure_depth * .25);
138  const G4double edge2_half_depth = sqrt(width_adj2 * width_adj2 + enclosure_depth * enclosure_depth * .25);
139 
140  // projective center
141  const G4double half_projection_ratio = 0.5 * (-width_adj1 + width_adj2) / enclosure_half_height_half_width;
142  const G4double projection_center_y = enclosure_center - ((enclosure_depth * 0.5) / half_projection_ratio);
143  const G4double projection_center_x = center_adj / half_projection_ratio;
144 
145  // blocks azimuthal segmentation
146  const int phi_bin_in_sec = get_geom_v3()->get_max_phi_bin_in_sec();
147  assert(phi_bin_in_sec >= 1);
148  const G4double block_azimuth_angle = (edge2_tilt_angle - edge1_tilt_angle) / phi_bin_in_sec;
149  assert(block_azimuth_angle > 0);
150  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()))
151  {
152  std::cout << "angle/nsec out of range: " << M_PI * std::numeric_limits<G4double>::epsilon() << std::endl;
153  gSystem->Exit(1);
154  }
155  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);
156  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);
157  G4double block_width_ratio = 0;
158  for (int sa = 0; sa < phi_bin_in_sec; ++sa)
159  {
160  block_width_ratio += 1 / cos(block_azimuth_angle * (0.5 + sa) + edge1_tilt_angle);
161  }
162  const G4double block_half_height_width = (block_edge1_half_width + block_edge2_half_width) / block_width_ratio;
163  assert(block_half_height_width > 0);
164 
165  // write out the azimuthal block geometry
166  // block azimuth geometry records
167  struct block_azimuth_geom
168  {
169  G4double angle;
170  G4double projection_center_y;
171  G4double projection_center_x;
172  G4double projection_length;
173  };
174  std::vector<block_azimuth_geom> block_azimuth_geoms(phi_bin_in_sec,
175  block_azimuth_geom{
176  std::numeric_limits<double>::quiet_NaN(),
177  std::numeric_limits<double>::quiet_NaN(),
178  std::numeric_limits<double>::quiet_NaN(),
179  std::numeric_limits<double>::quiet_NaN()}); // [phi-bin in sector] -> azimuth geometry
180  G4double block_x_edge1 = block_edge1_half_width;
181  for (int sa = 0; sa < phi_bin_in_sec; ++sa)
182  {
183  block_azimuth_geom& geom = block_azimuth_geoms[sa];
184 
185  geom.angle = block_azimuth_angle * (0.5 + sa) + edge1_tilt_angle;
186  const G4double block_x_size = block_half_height_width / cos(geom.angle);
187  assert(block_x_size > 0);
188  const G4double x_center = block_x_edge1 - 0.5 * block_x_size;
189 
190  // projection center per block
191  geom.projection_length = block_half_height_width / 2. / tan(block_azimuth_angle / 2.);
192  assert(geom.projection_length > 0);
193  geom.projection_center_y = enclosure_center - geom.projection_length * cos(geom.angle);
194  geom.projection_center_x = x_center + geom.projection_length * sin(geom.angle);
195 
196  // next step
197  block_x_edge1 -= block_x_size;
198  }
199 
200  // write out the azimuthal block divider's geometry
201  struct block_divider_azimuth_geom
202  {
203  G4double angle;
204  G4double projection_center_y;
205  G4double projection_center_x;
206  G4double thickness; // thickness in the approximate azimuth direction
207  G4double radial_displacement;
208  G4double width;
209  };
210  assert(phi_bin_in_sec >= 1);
211  std::vector<block_divider_azimuth_geom> divider_azimuth_geoms(phi_bin_in_sec - 1,
212  block_divider_azimuth_geom{
213  std::numeric_limits<double>::quiet_NaN(),
214  std::numeric_limits<double>::quiet_NaN(),
215  std::numeric_limits<double>::quiet_NaN(),
216  std::numeric_limits<double>::quiet_NaN(),
217  std::numeric_limits<double>::quiet_NaN(),
218  std::numeric_limits<double>::quiet_NaN()});
219 
220  if (get_geom_v3()->get_sidewall_thickness() > 0)
221  {
222  for (int sa = 0; sa < phi_bin_in_sec - 1; ++sa)
223  {
224  block_divider_azimuth_geom& geom = divider_azimuth_geoms[sa];
225 
226  geom.angle = 0.5 * (block_azimuth_geoms[sa].angle + block_azimuth_geoms[sa + 1].angle);
227  geom.projection_center_y = 0.5 * (block_azimuth_geoms[sa].projection_center_y + block_azimuth_geoms[sa + 1].projection_center_y);
228  geom.projection_center_x = 0.5 * (block_azimuth_geoms[sa].projection_center_x + block_azimuth_geoms[sa + 1].projection_center_x);
229  geom.radial_displacement = 0.5 * (block_azimuth_geoms[sa].projection_length + block_azimuth_geoms[sa + 1].projection_length);
230 
231  geom.thickness = 2.0 * get_geom_v3()->get_assembly_spacing() * cm * cos(block_azimuth_angle / 2.) - 2 * um;
232  geom.width = get_geom_v3()->get_divider_width() * cm;
233  }
234  }
235 
236  if (fabs(block_x_edge1 - (-block_edge2_half_width)) > get_geom_v3()->get_assembly_spacing() * cm)
237  {
238  std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg - ERROR - " << std::endl
239  << "\t block_x_edge1 = " << block_x_edge1 << std::endl
240  << "\t block_edge2_half_width = " << block_edge2_half_width << std::endl
241  << "\t fabs(block_x_edge1 - (-block_edge2_half_width)) = " << fabs(block_x_edge1 - (-block_edge2_half_width)) << std::endl
242  << "\t get_geom_v3()->get_assembly_spacing() * cm = " << get_geom_v3()->get_assembly_spacing() * cm << std::endl;
243  }
244  if (!(fabs(block_x_edge1 - (-block_edge2_half_width)) < get_geom_v3()->get_assembly_spacing() * cm)) // closure check
245  {
246  std::cout << "closure check failed: " << fabs(block_x_edge1 - (-block_edge2_half_width)) << std::endl;
247  gSystem->Exit(1);
248  }
249 
250  if (Verbosity())
251  {
252  std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg - " << std::endl
253  << "\t edge1_tilt_angle = " << edge1_tilt_angle << std::endl
254  << "\t edge2_tilt_angle = " << edge2_tilt_angle << std::endl
255  << "\t projection_center_y = " << projection_center_y << std::endl
256  << "\t projection_center_x = " << projection_center_x << std::endl
257  << "\t block_azimuth_angle = " << block_azimuth_angle << std::endl
258  << "\t block_edge1_half_width = " << block_edge1_half_width << std::endl
259  << "\t block_edge2_half_width = " << block_edge2_half_width << std::endl
260  << "\t block_width_ratio = " << block_width_ratio << std::endl
261  << "\t block_half_height_width = " << block_half_height_width << std::endl;
262 
263  for (int sa = 0; sa < phi_bin_in_sec; ++sa)
264  {
265  std::cout << "\t block[" << sa << "].angle = " << block_azimuth_geoms[sa].angle << std::endl;
266  std::cout << "\t block[" << sa << "].projection_center_y = " << block_azimuth_geoms[sa].projection_center_y << std::endl;
267  std::cout << "\t block[" << sa << "].projection_center_x = " << block_azimuth_geoms[sa].projection_center_x << std::endl;
268  }
269  for (int sa = 0; sa < phi_bin_in_sec - 1; ++sa)
270  {
271  std::cout << "\t divider[" << sa << "].angle = " << divider_azimuth_geoms[sa].angle << std::endl;
272  std::cout << "\t divider[" << sa << "].projection_center_x = " << divider_azimuth_geoms[sa].projection_center_x << std::endl;
273  std::cout << "\t divider[" << sa << "].projection_center_y = " << divider_azimuth_geoms[sa].projection_center_y << std::endl;
274  std::cout << "\t divider[" << sa << "].radial_displacement = " << divider_azimuth_geoms[sa].radial_displacement << std::endl;
275  std::cout << "\t divider[" << sa << "].thickness = " << divider_azimuth_geoms[sa].thickness << std::endl;
276  std::cout << "\t divider[" << sa << "].width = " << divider_azimuth_geoms[sa].width << std::endl;
277  }
278  }
279 
280  assert(enclosure_depth > 10 * cm);
281 
282  G4VSolid* sec_solid = new G4Trap(
283  /*const G4String& pName*/ G4String(GetName() + std::string("_sec_trap")),
284  enclosure_depth * 0.5, // G4double pDz,
285  center_tilt_angle, halfpi, // G4double pTheta, G4double pPhi,
286  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,
287  0, // G4double pAlp1,
288  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,
289  0 // G4double pAlp2 //
290  );
291  G4Transform3D sec_solid_transform = G4TranslateY3D(enclosure_center) * G4RotateY3D(halfpi) * G4RotateX3D(-halfpi);
292  G4VSolid* sec_solid_place = new G4DisplacedSolid(G4String(GetName() + std::string("_sec")), sec_solid, sec_solid_transform);
293 
294  recoConsts* rc = recoConsts::instance();
295  G4Material* cylinder_mat = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
296  assert(cylinder_mat);
297 
298  G4LogicalVolume* sec_logic = new G4LogicalVolume(sec_solid_place, cylinder_mat,
299  G4String(G4String(GetName() + std::string("_sec"))), nullptr, nullptr, nullptr);
300 
301  GetDisplayAction()->AddVolume(sec_logic, "Sector");
302 
303  // construct walls
304 
305  G4Material* wall_mat = GetDetectorMaterial(get_geom_v3()->get_sidewall_mat());
306  assert(wall_mat);
307 
308  if (get_geom_v3()->get_sidewall_thickness() > 0)
309  {
310  // end walls
311  if (get_geom_v3()->get_construction_verbose() >= 1)
312  {
313  std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg::" << GetName()
314  << " - construct end walls." << std::endl;
315  }
316  // G4Tubs* wall_solid = new G4Tubs(G4String(GetName() + std::string("_EndWall")),
317  // get_geom_v3()->get_radius() * cm + get_geom_v3()->get_sidewall_outer_torr() * cm,
318  // get_geom_v3()->get_max_radius() * cm - get_geom_v3()->get_sidewall_outer_torr() * cm,
319  // get_geom_v3()->get_sidewall_thickness() * cm / 2.0,
320  // halfpi - pi / get_geom_v3()->get_azimuthal_n_sec(),
321  // twopi / get_geom_v3()->get_azimuthal_n_sec());
322  const G4double side_wall_half_thickness = get_geom_v3()->get_sidewall_thickness() * cm / 2.0;
323  G4VSolid* wall_solid = new G4Trap(G4String(GetName() + std::string("_EndWall_trap")),
324  enclosure_depth * 0.5, // G4double pDz,
325  center_tilt_angle, halfpi, // G4double pTheta, G4double pPhi,
326  inner_half_width, side_wall_half_thickness, side_wall_half_thickness, // G4double pDy1, G4double pDx1, G4double pDx2,
327  0, // G4double pAlp1,
328  outter_half_width, side_wall_half_thickness, side_wall_half_thickness, // G4double pDy2, G4double pDx3, G4double pDx4,
329  0 // G4double pAlp2 //
330  );
331  G4VSolid* wall_solid_place = new G4DisplacedSolid(G4String(GetName() + std::string("_EndWall")), wall_solid, sec_solid_transform);
332 
333  G4LogicalVolume* wall_logic = new G4LogicalVolume(wall_solid_place, wall_mat,
334  G4String(G4String(GetName() + std::string("_EndWall"))), nullptr, nullptr,
335  nullptr);
336  GetDisplayAction()->AddVolume(wall_logic, "Wall");
337 
338  using z_locations_t = std::map<int, double>;
339  z_locations_t z_locations;
340  z_locations[000] = get_geom_v3()->get_sidewall_thickness() * cm / 2.0 + get_geom_v3()->get_assembly_spacing() * cm;
341  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);
342  z_locations[100] = -(get_geom_v3()->get_sidewall_thickness() * cm / 2.0 + get_geom_v3()->get_assembly_spacing() * cm);
343  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));
344 
345  for (z_locations_t::value_type& val : z_locations)
346  {
347  if (get_geom_v3()->get_construction_verbose() >= 2)
348  {
349  std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg::"
350  << GetName() << " - constructed End Wall ID " << val.first
351  << " @ Z = " << val.second << std::endl;
352  }
353  G4Transform3D wall_trans = G4TranslateZ3D(val.second);
354 
355  G4PVPlacement* wall_phys = new G4PVPlacement(wall_trans, wall_logic,
356  G4String(GetName()) + G4String("_EndWall_") + std::to_string(val.first), sec_logic,
357  false, val.first, OverlapCheck());
358 
359  calo_vol[wall_phys] = val.first;
360  assert(gdml_config);
361  gdml_config->exclude_physical_vol(wall_phys);
362  }
363  }
364  //
365  if (get_geom_v3()->get_sidewall_thickness() > 0)
366  {
367  // side walls
368  if (get_geom_v3()->get_construction_verbose() >= 1)
369  {
370  std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg::" << GetName()
371  << " - construct side walls." << std::endl;
372  }
373 
374  using sign_t = std::map<int, std::pair<int, int>>;
375  sign_t signs;
376  signs[100] = std::make_pair(+1, +1);
377  signs[101] = std::make_pair(+1, -1);
378  signs[200] = std::make_pair(-1, +1);
379  signs[201] = std::make_pair(-1, -1);
380 
381  for (sign_t::value_type& val : signs)
382  {
383  const int sign_z = val.second.first;
384  const int sign_azimuth = val.second.second;
385 
386  if (get_geom_v3()->get_construction_verbose() >= 2)
387  {
388  std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg::"
389  << GetName() << " - constructed Side Wall ID " << val.first
390  << " with"
391  << " Shift X = "
392  << sign_azimuth * (get_geom_v3()->get_sidewall_thickness() * cm / 2.0 + get_geom_v3()->get_sidewall_outer_torr() * cm)
393  << " Rotation Z = "
394  << sign_azimuth * pi / get_geom_v3()->get_azimuthal_n_sec()
395  << " Shift Z = " << sign_z * (get_geom_v3()->get_length() * cm / 4)
396  << std::endl;
397  }
398  const G4double azimuth_roate = sign_azimuth > 0 ? edge1_tilt_angle : edge2_tilt_angle;
399  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);
400 
401  G4Box* wall_solid = new G4Box(G4String(GetName() + G4String("_SideWall_") + std::to_string(val.first)),
402  get_geom_v3()->get_sidewall_thickness() * cm / 2.0,
403  edge_half_depth,
404  (get_geom_v3()->get_length() / 2. - 2 * (get_geom_v3()->get_sidewall_thickness() + 2. * get_geom_v3()->get_assembly_spacing())) * cm * .5);
405 
406  G4LogicalVolume* wall_logic = new G4LogicalVolume(wall_solid, wall_mat,
407  G4String(G4String(GetName() + G4String("_SideWall_") + std::to_string(val.first))), nullptr, nullptr,
408  nullptr);
409  GetDisplayAction()->AddVolume(wall_logic, "Wall");
410 
411  const G4Transform3D wall_trans =
412  G4TranslateZ3D(sign_z * (get_geom_v3()->get_length() * cm / 4)) *
413  G4TranslateY3D(enclosure_center) *
414  G4TranslateX3D(sign_azimuth * enclosure_half_height_half_width) *
415  G4RotateZ3D(azimuth_roate) *
416  G4TranslateX3D(-sign_azimuth * (get_geom_v3()->get_sidewall_thickness() * cm / 2.0 + get_geom_v3()->get_sidewall_outer_torr() * cm));
417 
418  G4PVPlacement* wall_phys = new G4PVPlacement(wall_trans, wall_logic,
419  G4String(GetName()) + G4String("_SideWall_") + std::to_string(val.first), sec_logic,
420  false, val.first, OverlapCheck());
421 
422  calo_vol[wall_phys] = val.first;
423  assert(gdml_config);
424  gdml_config->exclude_physical_vol(wall_phys);
425  }
426  }
427 
428  // construct dividers
429  if (get_geom_v3()->get_divider_width() > 0)
430  {
431  if (get_geom_v3()->get_construction_verbose() >= 1)
432  {
433  std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg::" << GetName()
434  << " - construct dividers" << std::endl;
435  }
436 
437  G4Material* divider_mat = GetDetectorMaterial(get_geom_v3()->get_divider_mat());
438  assert(divider_mat);
439 
440  int ID = 300;
441  for (const auto& geom : divider_azimuth_geoms)
442  {
443  G4Box* divider_solid = new G4Box(G4String(GetName() + G4String("_Divider_") + std::to_string(ID)),
444  geom.thickness / 2.0,
445  geom.width / 2.,
446  (get_geom_v3()->get_length() / 2. - 2 * (get_geom_v3()->get_sidewall_thickness() + 2. * get_geom_v3()->get_assembly_spacing())) * cm * .5);
447 
448  G4LogicalVolume* wall_logic = new G4LogicalVolume(divider_solid, divider_mat,
449  G4String(G4String(GetName() + G4String("_Divider_") + std::to_string(ID))), nullptr, nullptr,
450  nullptr);
451  GetDisplayAction()->AddVolume(wall_logic, "Divider");
452 
453  for (int sign_z = -1; sign_z <= 1; sign_z += 2)
454  {
455  G4Transform3D wall_trans =
456  G4TranslateX3D(geom.projection_center_x) *
457  G4TranslateY3D(geom.projection_center_y) *
458  G4RotateZ3D(geom.angle) *
459  G4TranslateY3D(geom.radial_displacement) *
460  G4TranslateZ3D(sign_z * (get_geom_v3()->get_length() * cm / 4));
461 
462  G4PVPlacement* wall_phys = new G4PVPlacement(wall_trans, wall_logic,
463  G4String(GetName()) + G4String("_Divider_") + std::to_string(ID), sec_logic,
464  false, ID, OverlapCheck());
465 
466  calo_vol[wall_phys] = ID;
467  assert(gdml_config);
468  gdml_config->exclude_physical_vol(wall_phys);
469 
470  if (Verbosity())
471  {
472  std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg - placing divider " << wall_phys->GetName() << " copy ID " << ID << std::endl;
473  }
474 
475  ++ID;
476  }
477  } // for (const auto & geom : divider_azimuth_geoms)
478  }
479 
480  // // construct towers
481  //
482  for (const SpacalGeom_t::tower_map_t::value_type& val : get_geom_v3()->get_sector_tower_map())
483  {
484  SpacalGeom_t::geom_tower g_tower = val.second;
485 
486  const int tower_id = g_tower.id;
487  const int tower_phi_id_in_sec = tower_id % 10;
488  assert(tower_phi_id_in_sec >= 0);
489  assert(tower_phi_id_in_sec < phi_bin_in_sec);
490 
491  const auto& block_azimuth_geom = block_azimuth_geoms.at(tower_phi_id_in_sec);
492 
493  G4LogicalVolume* LV_tower = Construct_Tower(g_tower);
494 
495  G4Transform3D block_trans =
496  G4TranslateX3D(block_azimuth_geom.projection_center_x) *
497  G4TranslateY3D(block_azimuth_geom.projection_center_y) *
498  G4RotateZ3D(block_azimuth_geom.angle) *
499  G4TranslateX3D(g_tower.centralX * cm) *
500  G4TranslateY3D(g_tower.centralY * cm) *
501  G4TranslateZ3D(g_tower.centralZ * cm) *
502  G4RotateX3D(g_tower.pRotationAngleX * rad);
503 
504  const bool overlapcheck_block = OverlapCheck() and (get_geom_v3()->get_construction_verbose() >= 2);
505 
506  G4PVPlacement* block_phys = new G4PVPlacement(block_trans, LV_tower,
507  G4String(GetName()) + G4String("_Tower_") + std::to_string(g_tower.id), sec_logic, false,
508  g_tower.id, overlapcheck_block);
509  block_vol[block_phys] = g_tower.id;
510  assert(gdml_config);
511  gdml_config->exclude_physical_vol(block_phys);
512 
513  if (g_tower.LightguideHeight > 0)
514  {
515  // also build a light guide
516 
517  for (int ix = 0; ix < g_tower.NSubtowerX; ix++)
518  // int ix = 0;
519  {
520  for (int iy = 0; iy < g_tower.NSubtowerY; iy++)
521  // int iy = 0;
522  {
523  G4LogicalVolume* LV_lg = Construct_LightGuide(g_tower, ix,
524  iy);
525 
526  G4PVPlacement* lg_phys = new G4PVPlacement(block_trans, LV_lg, LV_lg->GetName(),
527  sec_logic, false, g_tower.id, overlapcheck_block);
528 
529  block_vol[lg_phys] = g_tower.id * 100 + ix * 10 + iy;
530 
531  assert(gdml_config);
532  gdml_config->exclude_physical_vol(lg_phys);
533  }
534  }
535  }
536  }
537 
538  std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_AzimuthalSeg::" << GetName()
539  << " - constructed " << get_geom_v3()->get_sector_tower_map().size()
540  << " unique towers" << std::endl;
541 
542  return std::make_pair(sec_logic, G4Transform3D::Identity);
543 }
544 
546 int PHG4FullProjTiltedSpacalDetector::Construct_Fibers_SameLengthFiberPerTower(
547  const PHG4FullProjTiltedSpacalDetector::SpacalGeom_t::geom_tower& g_tower,
548  G4LogicalVolume* LV_tower)
549 {
550  // construct fibers
551 
552  // first check out the fibers geometry
553 
554  using fiber_par_map = std::map<int, std::pair<G4Vector3D, G4Vector3D>>;
555  fiber_par_map fiber_par;
556  G4double min_fiber_length = g_tower.pDz * cm * 4;
557 
558  G4Vector3D v_zshift = G4Vector3D(tan(g_tower.pTheta) * cos(g_tower.pPhi),
559  tan(g_tower.pTheta) * sin(g_tower.pPhi), 1) *
560  g_tower.pDz;
561  // int fiber_ID = 0;
562  for (int ix = 0; ix < g_tower.NFiberX; ix++)
563  // int ix = 0;
564  {
565  const double weighted_ix = static_cast<double>(ix) / (g_tower.NFiberX - 1.);
566 
567  const double weighted_pDx1 = (g_tower.pDx1 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
568  const double weighted_pDx2 = (g_tower.pDx2 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
569 
570  const double weighted_pDx3 = (g_tower.pDx3 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
571  const double weighted_pDx4 = (g_tower.pDx4 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
572 
573  for (int iy = 0; iy < g_tower.NFiberY; iy++)
574  // int iy = 0;
575  {
576  if ((ix + iy) % 2 == 1)
577  continue; // make a triangle pattern
578 
579  const double weighted_iy = static_cast<double>(iy) / (g_tower.NFiberY - 1.);
580 
581  const double weighted_pDy1 = (g_tower.pDy1 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
582  const double weighted_pDy2 = (g_tower.pDy2 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
583 
584  const double weighted_pDx12 = weighted_pDx1 * (1 - weighted_iy) + weighted_pDx2 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp1);
585  const double weighted_pDx34 = weighted_pDx3 * (1 - weighted_iy) + weighted_pDx4 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp2);
586 
587  G4Vector3D v1 = G4Vector3D(weighted_pDx12, weighted_pDy1, 0) - v_zshift;
588  G4Vector3D v2 = G4Vector3D(weighted_pDx34, weighted_pDy2, 0) + v_zshift;
589 
590  G4Vector3D vector_fiber = (v2 - v1);
591  vector_fiber *= (vector_fiber.mag() - get_geom_v3()->get_fiber_outer_r()) / vector_fiber.mag(); // shrink by fiber boundary protection
592  G4Vector3D center_fiber = (v2 + v1) / 2;
593 
594  // convert to Geant4 units
595  vector_fiber *= cm;
596  center_fiber *= cm;
597 
598  const int fiber_ID = g_tower.compose_fiber_id(ix, iy);
599  fiber_par[fiber_ID] = std::make_pair(vector_fiber,
600  center_fiber);
601 
602  const G4double fiber_length = vector_fiber.mag();
603 
604  min_fiber_length = std::min(fiber_length, min_fiber_length);
605 
606  // ++fiber_ID;
607  }
608  }
609 
610  int fiber_count = 0;
611 
612  const G4double fiber_length = min_fiber_length;
613  std::vector<G4double> fiber_cut;
614 
615  std::stringstream ss;
616  ss << std::string("_Tower") << g_tower.id;
617  G4LogicalVolume* fiber_logic = Construct_Fiber(fiber_length, ss.str());
618 
619  for (const fiber_par_map::value_type& val : fiber_par)
620  {
621  const int fiber_ID = val.first;
622  G4Vector3D vector_fiber = val.second.first;
623  G4Vector3D center_fiber = val.second.second;
624  const G4double optimal_fiber_length = vector_fiber.mag();
625 
626  const G4Vector3D v1 = center_fiber - 0.5 * vector_fiber;
627 
628  // keep a statistics
629  assert(optimal_fiber_length - fiber_length >= 0);
630  fiber_cut.push_back(optimal_fiber_length - fiber_length);
631 
632  center_fiber += (fiber_length / optimal_fiber_length - 1) * 0.5 * vector_fiber;
633  vector_fiber *= fiber_length / optimal_fiber_length;
634 
635  // const G4Vector3D v1_new = center_fiber - 0.5 *vector_fiber;
636 
637  if (get_geom_v3()->get_construction_verbose() >= 3)
638  std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_Fibers_SameLengthFiberPerTower::" << GetName()
639  << " - constructed fiber " << fiber_ID << ss.str() //
640  << ", Length = " << optimal_fiber_length << "-"
641  << (optimal_fiber_length - fiber_length) << "mm, " //
642  << "x = " << center_fiber.x() << "mm, " //
643  << "y = " << center_fiber.y() << "mm, " //
644  << "z = " << center_fiber.z() << "mm, " //
645  << "vx = " << vector_fiber.x() << "mm, " //
646  << "vy = " << vector_fiber.y() << "mm, " //
647  << "vz = " << vector_fiber.z() << "mm, " //
648  << std::endl;
649 
650  const G4double rotation_angle = G4Vector3D(0, 0, 1).angle(vector_fiber);
651  const G4Vector3D rotation_axis =
652  rotation_angle == 0 ? G4Vector3D(1, 0, 0) : G4Vector3D(0, 0, 1).cross(vector_fiber);
653 
654  G4Transform3D fiber_place(
655  G4Translate3D(center_fiber.x(), center_fiber.y(), center_fiber.z()) * G4Rotate3D(rotation_angle, rotation_axis));
656 
657  std::stringstream name;
658  name << GetName() + std::string("_Tower") << g_tower.id << "_fiber"
659  << ss.str();
660 
661  const bool overlapcheck_fiber = OverlapCheck() and (get_geom_v3()->get_construction_verbose() >= 3);
662  G4PVPlacement* fiber_physi = new G4PVPlacement(fiber_place, fiber_logic,
663  G4String(name.str()), LV_tower, false, fiber_ID,
664  overlapcheck_fiber);
665  fiber_vol[fiber_physi] = fiber_ID;
666  assert(gdml_config);
667  gdml_config->exclude_physical_vol(fiber_physi);
668 
669  fiber_count++;
670  }
671 
672  if (get_geom_v3()->get_construction_verbose() >= 2)
673  std::cout
674  << "PHG4FullProjTiltedSpacalDetector::Construct_Fibers_SameLengthFiberPerTower::"
675  << GetName() << " - constructed tower ID " << g_tower.id << " with "
676  << fiber_count << " fibers. Average fiber length cut = "
677  << std::accumulate(fiber_cut.begin(), fiber_cut.end(), 0.0) / fiber_cut.size() << " mm" << std::endl;
678 
679  return fiber_count;
680 }
681 
683 int PHG4FullProjTiltedSpacalDetector::Construct_Fibers(
684  const PHG4FullProjTiltedSpacalDetector::SpacalGeom_t::geom_tower& g_tower,
685  G4LogicalVolume* LV_tower)
686 {
687  G4Vector3D v_zshift = G4Vector3D(tan(g_tower.pTheta) * cos(g_tower.pPhi),
688  tan(g_tower.pTheta) * sin(g_tower.pPhi), 1) *
689  g_tower.pDz;
690  int fiber_cnt = 0;
691  for (int ix = 0; ix < g_tower.NFiberX; ix++)
692  {
693  const double weighted_ix = static_cast<double>(ix) / (g_tower.NFiberX - 1.);
694 
695  const double weighted_pDx1 = (g_tower.pDx1 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
696  const double weighted_pDx2 = (g_tower.pDx2 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
697 
698  const double weighted_pDx3 = (g_tower.pDx3 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
699  const double weighted_pDx4 = (g_tower.pDx4 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
700 
701  for (int iy = 0; iy < g_tower.NFiberY; iy++)
702  {
703  if ((ix + iy) % 2 == 1)
704  continue; // make a triangle pattern
705  const int fiber_ID = g_tower.compose_fiber_id(ix, iy);
706 
707  const double weighted_iy = static_cast<double>(iy) / (g_tower.NFiberY - 1.);
708 
709  const double weighted_pDy1 = (g_tower.pDy1 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
710  const double weighted_pDy2 = (g_tower.pDy2 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
711 
712  const double weighted_pDx12 = weighted_pDx1 * (1 - weighted_iy) + weighted_pDx2 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp1);
713  const double weighted_pDx34 = weighted_pDx3 * (1 - weighted_iy) + weighted_pDx4 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp2);
714 
715  G4Vector3D v1 = G4Vector3D(weighted_pDx12, weighted_pDy1, 0) - v_zshift;
716  G4Vector3D v2 = G4Vector3D(weighted_pDx34, weighted_pDy2, 0) + v_zshift;
717 
718  G4Vector3D vector_fiber = (v2 - v1);
719  vector_fiber *= (vector_fiber.mag() - get_geom_v3()->get_fiber_outer_r()) / vector_fiber.mag(); // shrink by fiber boundary protection
720  G4Vector3D center_fiber = (v2 + v1) / 2;
721 
722  // convert to Geant4 units
723  vector_fiber *= cm;
724  center_fiber *= cm;
725 
726  const G4double fiber_length = vector_fiber.mag();
727 
728  std::stringstream ss;
729  ss << std::string("_Tower") << g_tower.id;
730  ss << "_x" << ix;
731  ss << "_y" << iy;
732  G4LogicalVolume* fiber_logic = Construct_Fiber(fiber_length,
733  ss.str());
734 
735  if (get_geom_v3()->get_construction_verbose() >= 3)
736  std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_Fibers::" << GetName()
737  << " - constructed fiber " << fiber_ID << ss.str() //
738  << ", Length = " << fiber_length << "mm, " //
739  << "x = " << center_fiber.x() << "mm, " //
740  << "y = " << center_fiber.y() << "mm, " //
741  << "z = " << center_fiber.z() << "mm, " //
742  << "vx = " << vector_fiber.x() << "mm, " //
743  << "vy = " << vector_fiber.y() << "mm, " //
744  << "vz = " << vector_fiber.z() << "mm, " //
745  << std::endl;
746 
747  const G4double rotation_angle = G4Vector3D(0, 0, 1).angle(
748  vector_fiber);
749  const G4Vector3D rotation_axis =
750  rotation_angle == 0 ? G4Vector3D(1, 0, 0) : G4Vector3D(0, 0, 1).cross(vector_fiber);
751 
752  G4Transform3D fiber_place(
753  G4Translate3D(center_fiber.x(), center_fiber.y(),
754  center_fiber.z()) *
755  G4Rotate3D(rotation_angle, rotation_axis));
756 
757  std::stringstream name;
758  name << GetName() + std::string("_Tower") << g_tower.id << "_fiber"
759  << ss.str();
760 
761  const bool overlapcheck_fiber = OverlapCheck() and (get_geom_v3()->get_construction_verbose() >= 3);
762  G4PVPlacement* fiber_physi = new G4PVPlacement(fiber_place,
763  fiber_logic, G4String(name.str()), LV_tower, false,
764  fiber_ID, overlapcheck_fiber);
765  fiber_vol[fiber_physi] = fiber_ID;
766  assert(gdml_config);
767  gdml_config->exclude_physical_vol(fiber_physi);
768 
769  ++fiber_cnt;
770  }
771  }
772 
773  if (get_geom_v3()->get_construction_verbose() >= 3)
774  std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_Fibers::" << GetName()
775  << " - constructed tower ID " << g_tower.id << " with " << fiber_cnt
776  << " fibers" << std::endl;
777 
778  return fiber_cnt;
779 }
780 
782 G4LogicalVolume*
783 PHG4FullProjTiltedSpacalDetector::Construct_Tower(
784  const PHG4FullProjTiltedSpacalDetector::SpacalGeom_t::geom_tower& g_tower)
785 {
786  std::stringstream sout;
787  sout << "_" << g_tower.id;
788  const G4String sTowerID(sout.str());
789 
790  // Processed PostionSeeds 1 from 1 1
791 
792  G4Trap* block_solid = new G4Trap(
793  /*const G4String& pName*/ G4String(GetName()) + sTowerID,
794  g_tower.pDz * cm, // G4double pDz,
795  g_tower.pTheta * rad, g_tower.pPhi * rad, // G4double pTheta, G4double pPhi,
796  g_tower.pDy1 * cm, g_tower.pDx1 * cm, g_tower.pDx2 * cm, // G4double pDy1, G4double pDx1, G4double pDx2,
797  g_tower.pAlp1 * rad, // G4double pAlp1,
798  g_tower.pDy2 * cm, g_tower.pDx3 * cm, g_tower.pDx4 * cm, // G4double pDy2, G4double pDx3, G4double pDx4,
799  g_tower.pAlp2 * rad // G4double pAlp2 //
800  );
801 
802  G4Material* cylinder_mat = GetDetectorMaterial(get_geom_v3()->get_absorber_mat());
803  assert(cylinder_mat);
804 
805  G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, cylinder_mat,
806  G4String(G4String(GetName()) + std::string("_Tower") + sTowerID), nullptr, nullptr,
807  nullptr);
808 
809  GetDisplayAction()->AddVolume(block_logic, "Block");
810 
811  // construct fibers
812 
813  if (get_geom_v3()->get_config() == SpacalGeom_t::kFullProjective_2DTaper_Tilted)
814  {
815  int fiber_count = Construct_Fibers(g_tower, block_logic);
816 
817  if (get_geom_v3()->get_construction_verbose() >= 2)
818  std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_Tower::" << GetName()
819  << " - constructed tower ID " << g_tower.id << " with "
820  << fiber_count << " fibers using Construct_Fibers" << std::endl;
821  }
822  else if (get_geom_v3()->get_config() == SpacalGeom_t::kFullProjective_2DTaper_Tilted_SameLengthFiberPerTower)
823  {
824  int fiber_count = Construct_Fibers_SameLengthFiberPerTower(g_tower,
825  block_logic);
826 
827  if (get_geom_v3()->get_construction_verbose() >= 2)
828  std::cout << "PHG4FullProjTiltedSpacalDetector::Construct_Tower::" << GetName()
829  << " - constructed tower ID " << g_tower.id << " with "
830  << fiber_count
831  << " fibers using Construct_Fibers_SameLengthFiberPerTower."
832  << "V = " << block_solid->GetCubicVolume() / (cm3) << "cm3, "
833  << "m = " << block_logic->GetMass() / gram << "gram, "
834  << "Density = " << (block_logic->GetMass() / gram) / (block_solid->GetCubicVolume() / cm3) << "g/cm3"
835  << std::endl;
836  }
837  else
838  {
839  G4ExceptionDescription message;
840  message << "can not recognize configuration type " << get_geom_v3()->get_config();
841 
842  G4Exception("PHG4FullProjTiltedSpacalDetector::Construct_Tower", "Wrong",
843  FatalException, message, "");
844  }
845 
846  return block_logic;
847 }
848 
849 G4LogicalVolume*
850 PHG4FullProjTiltedSpacalDetector::Construct_LightGuide(
851  const PHG4FullProjTiltedSpacalDetector::SpacalGeom_t::geom_tower& g_tower,
852  const int index_x, const int index_y)
853 {
854  assert(_geom);
855  std::stringstream sout;
856  sout << "_Lightguide_" << g_tower.id << "_" << index_x << "_" << index_y;
857  const G4String sTowerID(sout.str());
858 
859  assert(g_tower.LightguideHeight > 0);
860 
861  // light guide parameters in PHENIX units
862  const double weight_x1 = 1 - (double) index_y / g_tower.NSubtowerY;
863  const double weight_x2 = 1 - (double) (index_y + 1) / g_tower.NSubtowerY;
864  const double weight_xcenter = 1 - (double) (index_y + 0.5) / g_tower.NSubtowerY;
865 
866  assert(weight_x1 >= 0 and weight_x1 <= 1);
867  assert(weight_x2 >= 0 and weight_x2 <= 1);
868  assert(weight_xcenter >= 0 and weight_xcenter <= 1);
869 
870  const double lg_pDx1 = (g_tower.pDx1 * weight_x1 //
871  + g_tower.pDx2 * (1 - weight_x1)) /
872  g_tower.NSubtowerX;
873  const double lg_pDx2 = (g_tower.pDx1 * weight_x2 //
874  + g_tower.pDx2 * (1 - weight_x2)) /
875  g_tower.NSubtowerX;
876  const double lg_pDy1 = g_tower.pDy1 / g_tower.NSubtowerY;
877  const double lg_Alp1 = atan(
878  (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));
879 
880  const double shift_xcenter = (g_tower.pDx1 * weight_xcenter //
881  + g_tower.pDx2 * (1 - weight_xcenter)) //
882  * //
883  (-g_tower.NSubtowerX + 1. + 2 * index_x) / (double) (g_tower.NSubtowerX);
884  const double shift_ycenter = g_tower.pDy1 //
885  * //
886  (-g_tower.NSubtowerY + 1. + 2 * index_y) / (double) (g_tower.NSubtowerY);
887 
888  G4VSolid* block_solid = new G4Trap(
889  /*const G4String& pName*/ G4String(GetName()) + sTowerID,
890  0.5 * g_tower.LightguideHeight * cm, // G4double pDz,
891  0 * rad, 0 * rad, // G4double pTheta, G4double pPhi,
892  g_tower.LightguideTaperRatio * lg_pDy1 * cm,
893  g_tower.LightguideTaperRatio * lg_pDx1 * cm,
894  g_tower.LightguideTaperRatio * lg_pDx2 * cm, // G4double pDy1, G4double pDx1, G4double pDx2,
895  lg_Alp1 * rad, // G4double pAlp1,
896  lg_pDy1 * cm, lg_pDx1 * cm, lg_pDx2 * cm, // G4double pDy2, G4double pDx3, G4double pDx4,
897  lg_Alp1 * rad // G4double pAlp2 //
898  );
899 
900  block_solid = new G4DisplacedSolid(G4String(GetName() + "_displaced"),
901  block_solid, nullptr, //
902  G4ThreeVector( //
903  tan(g_tower.pTheta * rad) * cos(g_tower.pPhi * rad), //
904  tan(g_tower.pTheta * rad) * sin(g_tower.pPhi * rad), //
905  1) * // G4ThreeVector
906  -(g_tower.pDz) *
907  cm //
908  + G4ThreeVector(shift_xcenter * cm, shift_ycenter * cm, 0) // shit in subtower direction
909  + G4ThreeVector(0, 0, -0.5 * g_tower.LightguideHeight * cm) // shift in the light guide height
910  );
911 
912  G4Material* cylinder_mat = GetDetectorMaterial(g_tower.LightguideMaterial);
913  assert(cylinder_mat);
914 
915  G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, cylinder_mat,
916  G4String(G4String(GetName()) + std::string("_Tower") + sTowerID), nullptr, nullptr,
917  nullptr);
918 
919  GetDisplayAction()->AddMaterial("LightGuide", g_tower.LightguideMaterial);
920  GetDisplayAction()->AddVolume(block_logic, "LightGuide");
921 
922  return block_logic;
923 }
924 
925 void PHG4FullProjTiltedSpacalDetector::Print(const std::string& /*what*/) const
926 {
927  std::cout << "PHG4FullProjTiltedSpacalDetector::Print::" << GetName()
928  << " - Print Geometry:" << std::endl;
929  get_geom_v3()->Print();
930 
931  return;
932 }