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