Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4FullProjSpacalDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4FullProjSpacalDetector.cc
1 
9 #include "PHG4FullProjSpacalDetector.h"
10 
11 #include "PHG4SpacalDisplayAction.h"
12 
13 #include <g4gdml/PHG4GDMLConfig.hh>
14 
15 #include <phool/recoConsts.h>
16 
17 #include <Geant4/G4Box.hh>
18 #include <Geant4/G4Exception.hh> // for G4Exception, G4ExceptionD
19 #include <Geant4/G4ExceptionSeverity.hh> // for FatalException
20 #include <Geant4/G4LogicalVolume.hh>
21 #include <Geant4/G4PVPlacement.hh>
22 #include <Geant4/G4PhysicalConstants.hh>
23 #include <Geant4/G4String.hh> // for G4String
24 #include <Geant4/G4SystemOfUnits.hh>
25 #include <Geant4/G4Transform3D.hh>
26 #include <Geant4/G4Trap.hh>
27 #include <Geant4/G4Tubs.hh>
28 #include <Geant4/G4Types.hh> // for G4double
29 #include <Geant4/G4Vector3D.hh>
30 
31 #include <TSystem.h>
32 
33 #include <algorithm>
34 #include <cassert>
35 #include <cmath>
36 #include <iostream> // for operator<<, basic_ostream
37 #include <map> // for allocator, map<>::value_type
38 #include <numeric> // std::accumulate
39 #include <sstream>
40 #include <string> // std::string, std::to_string
41 #include <vector> // for vector
42 
43 class G4Material;
44 class PHCompositeNode;
45 
46 //_______________________________________________________________
47 // note this inactive thickness is ~1.5% of a radiation length
48 PHG4FullProjSpacalDetector::PHG4FullProjSpacalDetector(PHG4Subsystem* subsys, PHCompositeNode* Node,
49  const std::string& dnam, PHParameters* parameters, const int lyr)
50  : PHG4SpacalDetector(subsys, Node, dnam, parameters, lyr, false)
51 {
52  assert(_geom == nullptr);
53 
54  _geom = new SpacalGeom_t();
55  if (_geom == nullptr)
56  {
57  std::cout
58  << "PHG4FullProjSpacalDetector::Constructor - Fatal Error - invalid geometry object!"
59  << std::endl;
60  gSystem->Exit(1);
61  }
62 
63  // this class loads Chris Cullen 2D spacal design July 2015 by default.
64  // this step is deprecated now
65  // get_geom_v3()->load_demo_sector_tower_map_2015_Chris_Cullen_2D_spacal();
66 
67  assert(parameters);
68  get_geom_v3()->ImportParameters(*parameters);
69 
70  // std::cout <<"PHG4FullProjSpacalDetector::Constructor - get_geom_v3()->Print();"<<std::endl;
71  // get_geom_v3()->Print();
72 }
73 
74 //_______________________________________________________________
75 void PHG4FullProjSpacalDetector::ConstructMe(G4LogicalVolume* logicWorld)
76 {
77  if (get_geom_v3()->get_construction_verbose() >= 1)
78  {
79  std::cout << "PHG4FullProjSpacalDetector::Construct::" << GetName()
80  << " - start with PHG4SpacalDetector::Construct()." << std::endl;
81  }
82 
83  PHG4SpacalDetector::ConstructMe(logicWorld);
84 
85  if (get_geom_v3()->get_construction_verbose() >= 1)
86  {
87  std::cout << "PHG4FullProjSpacalDetector::Construct::" << GetName()
88  << " - Completed." << std::endl;
89  }
90 }
91 
92 std::pair<G4LogicalVolume*, G4Transform3D>
93 PHG4FullProjSpacalDetector::Construct_AzimuthalSeg()
94 {
95  if (!(get_geom_v3()->get_azimuthal_n_sec() > 4))
96  {
97  std::cout << "azimuthal_n_sec <= 4: " << get_geom_v3()->get_azimuthal_n_sec() << std::endl;
98  gSystem->Exit(1);
99  }
100 
101  G4Tubs* sec_solid = new G4Tubs(G4String(GetName() + std::string("_sec")),
102  get_geom_v3()->get_radius() * cm, get_geom_v3()->get_max_radius() * cm,
103  get_geom_v3()->get_length() * cm / 2.0,
104  halfpi - pi / get_geom_v3()->get_azimuthal_n_sec(),
105  twopi / get_geom_v3()->get_azimuthal_n_sec());
106 
107  recoConsts* rc = recoConsts::instance();
108  G4Material* cylinder_mat = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
109  assert(cylinder_mat);
110 
111  G4LogicalVolume* sec_logic = new G4LogicalVolume(sec_solid, cylinder_mat,
112  G4String(G4String(GetName() + std::string("_sec"))), nullptr, nullptr);
113 
114  GetDisplayAction()->AddVolume(sec_logic, "Sector");
115 
116  // construct walls
117 
118  G4Material* wall_mat = GetDetectorMaterial(get_geom_v3()->get_sidewall_mat());
119  assert(wall_mat);
120 
121  if (get_geom_v3()->get_sidewall_thickness() > 0)
122  {
123  // end walls
124  if (get_geom_v3()->get_construction_verbose() >= 1)
125  {
126  std::cout << "PHG4FullProjSpacalDetector::Construct_AzimuthalSeg::" << GetName()
127  << " - construct end walls." << std::endl;
128  }
129  G4Tubs* wall_solid = new G4Tubs(G4String(GetName() + std::string("_EndWall")),
130  get_geom_v3()->get_radius() * cm + get_geom_v3()->get_sidewall_outer_torr() * cm,
131  get_geom_v3()->get_max_radius() * cm - get_geom_v3()->get_sidewall_outer_torr() * cm,
132  get_geom_v3()->get_sidewall_thickness() * cm / 2.0,
133  halfpi - pi / get_geom_v3()->get_azimuthal_n_sec(),
134  twopi / get_geom_v3()->get_azimuthal_n_sec());
135 
136  G4LogicalVolume* wall_logic = new G4LogicalVolume(wall_solid, wall_mat,
137  G4String(G4String(GetName() + std::string("_EndWall"))), nullptr, nullptr,
138  nullptr);
139  GetDisplayAction()->AddVolume(wall_logic, "WallProj");
140 
141  using z_locations_t = std::map<int, double>;
142  z_locations_t z_locations;
143  z_locations[1000] = get_geom_v3()->get_sidewall_thickness() * cm / 2.0 + get_geom_v3()->get_assembly_spacing() * cm;
144  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);
145  z_locations[1100] = -(get_geom_v3()->get_sidewall_thickness() * cm / 2.0 + get_geom_v3()->get_assembly_spacing() * cm);
146  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));
147 
148  for (z_locations_t::value_type& val : z_locations)
149  {
150  if (get_geom_v3()->get_construction_verbose() >= 2)
151  std::cout << "PHG4FullProjSpacalDetector::Construct_AzimuthalSeg::"
152  << GetName() << " - constructed End Wall ID " << val.first
153  << " @ Z = " << val.second << std::endl;
154 
155  G4Transform3D wall_trans = G4TranslateZ3D(val.second);
156 
157  G4PVPlacement* wall_phys = new G4PVPlacement(wall_trans, wall_logic,
158  G4String(GetName()) + G4String("_EndWall"), sec_logic,
159  false, val.first, OverlapCheck());
160 
161  calo_vol[wall_phys] = val.first;
162  assert(gdml_config);
163  gdml_config->exclude_physical_vol(wall_phys);
164  }
165  }
166 
167  if (get_geom_v3()->get_sidewall_thickness() > 0)
168  {
169  // side walls
170  if (get_geom_v3()->get_construction_verbose() >= 1)
171  {
172  std::cout << "PHG4FullProjSpacalDetector::Construct_AzimuthalSeg::" << GetName()
173  << " - construct side walls." << std::endl;
174  }
175  G4Box* wall_solid = new G4Box(G4String(GetName() + std::string("_SideWall")),
176  get_geom_v3()->get_sidewall_thickness() * cm / 2.0,
177  get_geom_v3()->get_thickness() * cm / 2. - 2 * get_geom_v3()->get_sidewall_outer_torr() * cm,
178  (get_geom_v3()->get_length() / 2. - 2 * (get_geom_v3()->get_sidewall_thickness() + 2. * get_geom_v3()->get_assembly_spacing())) * cm * .5);
179 
180  G4LogicalVolume* wall_logic = new G4LogicalVolume(wall_solid, wall_mat,
181  G4String(G4String(GetName() + std::string("_SideWall"))), nullptr, nullptr,
182  nullptr);
183  GetDisplayAction()->AddVolume(wall_logic, "WallProj");
184 
185  using sign_t = std::map<int, std::pair<int, int>>;
186  sign_t signs;
187  signs[2000] = std::make_pair(+1, +1);
188  signs[2001] = std::make_pair(+1, -1);
189  signs[2100] = std::make_pair(-1, +1);
190  signs[2101] = std::make_pair(-1, -1);
191 
192  for (sign_t::value_type& val : signs)
193  {
194  const int sign_z = val.second.first;
195  const int sign_azimuth = val.second.second;
196 
197  if (get_geom_v3()->get_construction_verbose() >= 2)
198  std::cout << "PHG4FullProjSpacalDetector::Construct_AzimuthalSeg::"
199  << GetName() << " - constructed Side Wall ID " << val.first
200  << " with"
201  << " Shift X = "
202  << sign_azimuth * (get_geom_v3()->get_sidewall_thickness() * cm / 2.0 + get_geom_v3()->get_sidewall_outer_torr() * cm)
203  << " Rotation Z = "
204  << sign_azimuth * pi / get_geom_v3()->get_azimuthal_n_sec()
205  << " Shift Z = " << sign_z * (get_geom_v3()->get_length() * cm / 4)
206  << std::endl;
207 
208  G4Transform3D wall_trans = G4RotateZ3D(
209  sign_azimuth * pi / get_geom_v3()->get_azimuthal_n_sec()) *
210  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));
211 
212  G4PVPlacement* wall_phys = new G4PVPlacement(wall_trans, wall_logic,
213  G4String(GetName()) + G4String("_EndWall"), sec_logic,
214  false, val.first, OverlapCheck());
215 
216  calo_vol[wall_phys] = val.first;
217 
218  assert(gdml_config);
219  gdml_config->exclude_physical_vol(wall_phys);
220  }
221  }
222 
223  // construct towers
224 
225  for (const SpacalGeom_t::tower_map_t::value_type& val : get_geom_v3()->get_sector_tower_map())
226  {
227  const SpacalGeom_t::geom_tower& g_tower = val.second;
228  G4LogicalVolume* LV_tower = Construct_Tower(g_tower);
229 
230  G4Transform3D block_trans = G4TranslateX3D(g_tower.centralX * cm) * G4TranslateY3D(g_tower.centralY * cm) * G4TranslateZ3D(g_tower.centralZ * cm) * G4RotateX3D(g_tower.pRotationAngleX * rad);
231 
232  const bool overlapcheck_block = OverlapCheck() and (get_geom_v3()->get_construction_verbose() >= 2);
233 
234  G4PVPlacement* block_phys = new G4PVPlacement(block_trans, LV_tower,
235  G4String(GetName()) + G4String("_Tower"), sec_logic, false,
236  g_tower.id, overlapcheck_block);
237  block_vol[block_phys] = g_tower.id;
238 
239  assert(gdml_config);
240  gdml_config->exclude_physical_vol(block_phys);
241  }
242 
243  std::cout << "PHG4FullProjSpacalDetector::Construct_AzimuthalSeg::" << GetName()
244  << " - constructed " << get_geom_v3()->get_sector_tower_map().size()
245  << " unique towers" << std::endl;
246 
247  return std::make_pair(sec_logic, G4Transform3D::Identity);
248 }
249 
251 int PHG4FullProjSpacalDetector::Construct_Fibers_SameLengthFiberPerTower(
252  const PHG4FullProjSpacalDetector::SpacalGeom_t::geom_tower& g_tower,
253  G4LogicalVolume* LV_tower)
254 {
255  // construct fibers
256 
257  // first check out the fibers geometry
258 
259  using fiber_par_map = std::map<int, std::pair<G4Vector3D, G4Vector3D>>;
260  fiber_par_map fiber_par;
261  G4double min_fiber_length = g_tower.pDz * cm * 4;
262 
263  G4Vector3D v_zshift = G4Vector3D(tan(g_tower.pTheta) * cos(g_tower.pPhi),
264  tan(g_tower.pTheta) * sin(g_tower.pPhi), 1) *
265  g_tower.pDz;
266  // int fiber_ID = 0;
267  for (int ix = 0; ix < g_tower.NFiberX; ix++)
268  // int ix = 0;
269  {
270  const double weighted_ix = static_cast<double>(ix) / (g_tower.NFiberX - 1.);
271 
272  const double weighted_pDx1 = (g_tower.pDx1 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
273  const double weighted_pDx2 = (g_tower.pDx2 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
274 
275  const double weighted_pDx3 = (g_tower.pDx3 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
276  const double weighted_pDx4 = (g_tower.pDx4 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
277 
278  for (int iy = 0; iy < g_tower.NFiberY; iy++)
279  // int iy = 0;
280  {
281  if ((ix + iy) % 2 == 1)
282  continue; // make a triangle pattern
283 
284  const double weighted_iy = static_cast<double>(iy) / (g_tower.NFiberY - 1.);
285 
286  const double weighted_pDy1 = (g_tower.pDy1 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
287  const double weighted_pDy2 = (g_tower.pDy2 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
288 
289  const double weighted_pDx12 = weighted_pDx1 * (1 - weighted_iy) + weighted_pDx2 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp1);
290  const double weighted_pDx34 = weighted_pDx3 * (1 - weighted_iy) + weighted_pDx4 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp2);
291 
292  G4Vector3D v1 = G4Vector3D(weighted_pDx12, weighted_pDy1, 0) - v_zshift;
293  G4Vector3D v2 = G4Vector3D(weighted_pDx34, weighted_pDy2, 0) + v_zshift;
294 
295  G4Vector3D vector_fiber = (v2 - v1);
296  vector_fiber *= (vector_fiber.mag() - get_geom_v3()->get_fiber_outer_r()) / vector_fiber.mag(); // shrink by fiber boundary protection
297  G4Vector3D center_fiber = (v2 + v1) / 2;
298 
299  // convert to Geant4 units
300  vector_fiber *= cm;
301  center_fiber *= cm;
302 
303  const int fiber_ID = g_tower.compose_fiber_id(ix, iy);
304  fiber_par[fiber_ID] = std::make_pair(vector_fiber,
305  center_fiber);
306 
307  const G4double fiber_length = vector_fiber.mag();
308 
309  min_fiber_length = std::min(fiber_length, min_fiber_length);
310 
311  // ++fiber_ID;
312  }
313  }
314 
315  int fiber_count = 0;
316 
317  const G4double fiber_length = min_fiber_length;
318  std::vector<G4double> fiber_cut;
319 
320  std::stringstream ss;
321  ss << std::string("_Tower") << g_tower.id;
322  G4LogicalVolume* fiber_logic = Construct_Fiber(fiber_length, ss.str());
323 
324  for (const fiber_par_map::value_type& val : fiber_par)
325  {
326  const int fiber_ID = val.first;
327  G4Vector3D vector_fiber = val.second.first;
328  G4Vector3D center_fiber = val.second.second;
329  const G4double optimal_fiber_length = vector_fiber.mag();
330 
331  const G4Vector3D v1 = center_fiber - 0.5 * vector_fiber;
332 
333  // keep a statistics
334  assert(optimal_fiber_length - fiber_length >= 0);
335  fiber_cut.push_back(optimal_fiber_length - fiber_length);
336 
337  center_fiber += (fiber_length / optimal_fiber_length - 1) * 0.5 * vector_fiber;
338  vector_fiber *= fiber_length / optimal_fiber_length;
339 
340  // const G4Vector3D v1_new = center_fiber - 0.5 *vector_fiber;
341 
342  if (get_geom_v3()->get_construction_verbose() >= 3)
343  std::cout << "PHG4FullProjSpacalDetector::Construct_Fibers_SameLengthFiberPerTower::" << GetName()
344  << " - constructed fiber " << fiber_ID << ss.str() //
345  << ", Length = " << optimal_fiber_length << "-"
346  << (optimal_fiber_length - fiber_length) << "mm, " //
347  << "x = " << center_fiber.x() << "mm, " //
348  << "y = " << center_fiber.y() << "mm, " //
349  << "z = " << center_fiber.z() << "mm, " //
350  << "vx = " << vector_fiber.x() << "mm, " //
351  << "vy = " << vector_fiber.y() << "mm, " //
352  << "vz = " << vector_fiber.z() << "mm, " //
353  << std::endl;
354 
355  const G4double rotation_angle = G4Vector3D(0, 0, 1).angle(vector_fiber);
356  const G4Vector3D rotation_axis =
357  rotation_angle == 0 ? G4Vector3D(1, 0, 0) : G4Vector3D(0, 0, 1).cross(vector_fiber);
358 
359  G4Transform3D fiber_place(
360  G4Translate3D(center_fiber.x(), center_fiber.y(), center_fiber.z()) * G4Rotate3D(rotation_angle, rotation_axis));
361 
362  std::stringstream name;
363  name << GetName() + std::string("_Tower") << g_tower.id << "_fiber"
364  << ss.str();
365 
366  const bool overlapcheck_fiber = OverlapCheck() and (get_geom_v3()->get_construction_verbose() >= 3);
367  G4PVPlacement* fiber_physi = new G4PVPlacement(fiber_place, fiber_logic,
368  G4String(name.str()), LV_tower, false, fiber_ID,
369  overlapcheck_fiber);
370  fiber_vol[fiber_physi] = fiber_ID;
371 
372  assert(gdml_config);
373  gdml_config->exclude_physical_vol(fiber_physi);
374 
375  fiber_count++;
376  }
377 
378  if (get_geom_v3()->get_construction_verbose() >= 2)
379  std::cout
380  << "PHG4FullProjSpacalDetector::Construct_Fibers_SameLengthFiberPerTower::"
381  << GetName() << " - constructed tower ID " << g_tower.id << " with "
382  << fiber_count << " fibers. Average fiber length cut = "
383  << accumulate(fiber_cut.begin(), fiber_cut.end(), 0.0) / fiber_cut.size() << " mm" << std::endl;
384 
385  return fiber_count;
386 }
387 
389 int PHG4FullProjSpacalDetector::Construct_Fibers(
390  const PHG4FullProjSpacalDetector::SpacalGeom_t::geom_tower& g_tower,
391  G4LogicalVolume* LV_tower)
392 {
393  G4Vector3D v_zshift = G4Vector3D(tan(g_tower.pTheta) * cos(g_tower.pPhi),
394  tan(g_tower.pTheta) * sin(g_tower.pPhi), 1) *
395  g_tower.pDz;
396  int fiber_cnt = 0;
397  for (int ix = 0; ix < g_tower.NFiberX; ix++)
398  {
399  const double weighted_ix = static_cast<double>(ix) / (g_tower.NFiberX - 1.);
400 
401  const double weighted_pDx1 = (g_tower.pDx1 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
402  const double weighted_pDx2 = (g_tower.pDx2 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
403 
404  const double weighted_pDx3 = (g_tower.pDx3 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
405  const double weighted_pDx4 = (g_tower.pDx4 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
406 
407  for (int iy = 0; iy < g_tower.NFiberY; iy++)
408  {
409  if ((ix + iy) % 2 == 1)
410  continue; // make a triangle pattern
411  const int fiber_ID = g_tower.compose_fiber_id(ix, iy);
412 
413  const double weighted_iy = static_cast<double>(iy) / (g_tower.NFiberY - 1.);
414 
415  const double weighted_pDy1 = (g_tower.pDy1 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
416  const double weighted_pDy2 = (g_tower.pDy2 - g_tower.ModuleSkinThickness - get_geom_v3()->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
417 
418  const double weighted_pDx12 = weighted_pDx1 * (1 - weighted_iy) + weighted_pDx2 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp1);
419  const double weighted_pDx34 = weighted_pDx3 * (1 - weighted_iy) + weighted_pDx4 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp2);
420 
421  G4Vector3D v1 = G4Vector3D(weighted_pDx12, weighted_pDy1, 0) - v_zshift;
422  G4Vector3D v2 = G4Vector3D(weighted_pDx34, weighted_pDy2, 0) + v_zshift;
423 
424  G4Vector3D vector_fiber = (v2 - v1);
425  vector_fiber *= (vector_fiber.mag() - get_geom_v3()->get_fiber_outer_r()) / vector_fiber.mag(); // shrink by fiber boundary protection
426  G4Vector3D center_fiber = (v2 + v1) / 2;
427 
428  // convert to Geant4 units
429  vector_fiber *= cm;
430  center_fiber *= cm;
431 
432  const G4double fiber_length = vector_fiber.mag();
433 
434  std::stringstream ss;
435  ss << std::string("_Tower") << g_tower.id;
436  ss << "_x" << ix;
437  ss << "_y" << iy;
438  G4LogicalVolume* fiber_logic = Construct_Fiber(fiber_length,
439  ss.str());
440 
441  if (get_geom_v3()->get_construction_verbose() >= 3)
442  std::cout << "PHG4FullProjSpacalDetector::Construct_Fibers::" << GetName()
443  << " - constructed fiber " << fiber_ID << ss.str() //
444  << ", Length = " << fiber_length << "mm, " //
445  << "x = " << center_fiber.x() << "mm, " //
446  << "y = " << center_fiber.y() << "mm, " //
447  << "z = " << center_fiber.z() << "mm, " //
448  << "vx = " << vector_fiber.x() << "mm, " //
449  << "vy = " << vector_fiber.y() << "mm, " //
450  << "vz = " << vector_fiber.z() << "mm, " //
451  << std::endl;
452 
453  const G4double rotation_angle = G4Vector3D(0, 0, 1).angle(
454  vector_fiber);
455  const G4Vector3D rotation_axis =
456  rotation_angle == 0 ? G4Vector3D(1, 0, 0) : G4Vector3D(0, 0, 1).cross(vector_fiber);
457 
458  G4Transform3D fiber_place(
459  G4Translate3D(center_fiber.x(), center_fiber.y(),
460  center_fiber.z()) *
461  G4Rotate3D(rotation_angle, rotation_axis));
462 
463  std::stringstream name;
464  name << GetName() + std::string("_Tower") << g_tower.id << "_fiber"
465  << ss.str();
466 
467  const bool overlapcheck_fiber = OverlapCheck() and (get_geom_v3()->get_construction_verbose() >= 3);
468  G4PVPlacement* fiber_physi = new G4PVPlacement(fiber_place,
469  fiber_logic, G4String(name.str()), LV_tower, false,
470  fiber_ID, overlapcheck_fiber);
471  fiber_vol[fiber_physi] = fiber_ID;
472 
473  assert(gdml_config);
474  gdml_config->exclude_physical_vol(fiber_physi);
475 
476  ++fiber_cnt;
477  }
478  }
479 
480  if (get_geom_v3()->get_construction_verbose() >= 3)
481  std::cout << "PHG4FullProjSpacalDetector::Construct_Fibers::" << GetName()
482  << " - constructed tower ID " << g_tower.id << " with " << fiber_cnt
483  << " fibers" << std::endl;
484 
485  return fiber_cnt;
486 }
487 
489 G4LogicalVolume*
490 PHG4FullProjSpacalDetector::Construct_Tower(
491  const PHG4FullProjSpacalDetector::SpacalGeom_t::geom_tower& g_tower)
492 {
493  std::stringstream sout;
494  sout << "_" << g_tower.id;
495  const G4String sTowerID(sout.str());
496 
497  // Processed PostionSeeds 1 from 1 1
498 
499  G4Trap* block_solid = new G4Trap(
500  /*const G4String& pName*/ G4String(GetName()) + sTowerID,
501  g_tower.pDz * cm, // G4double pDz,
502  g_tower.pTheta * rad, g_tower.pPhi * rad, // G4double pTheta, G4double pPhi,
503  g_tower.pDy1 * cm, g_tower.pDx1 * cm, g_tower.pDx2 * cm, // G4double pDy1, G4double pDx1, G4double pDx2,
504  g_tower.pAlp1 * rad, // G4double pAlp1,
505  g_tower.pDy2 * cm, g_tower.pDx3 * cm, g_tower.pDx4 * cm, // G4double pDy2, G4double pDx3, G4double pDx4,
506  g_tower.pAlp2 * rad // G4double pAlp2 //
507  );
508 
509  G4Material* cylinder_mat = GetDetectorMaterial(get_geom_v3()->get_absorber_mat());
510  assert(cylinder_mat);
511 
512  G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, cylinder_mat,
513  G4String(G4String(GetName()) + std::string("_Tower") + sTowerID), nullptr, nullptr,
514  nullptr);
515 
516  GetDisplayAction()->AddVolume(block_logic, "Block");
517 
518  // construct fibers
519 
520  if (get_geom_v3()->get_config() == SpacalGeom_t::kFullProjective_2DTaper)
521  {
522  int fiber_count = Construct_Fibers(g_tower, block_logic);
523 
524  if (get_geom_v3()->get_construction_verbose() >= 2)
525  std::cout << "PHG4FullProjSpacalDetector::Construct_Tower::" << GetName()
526  << " - constructed tower ID " << g_tower.id << " with "
527  << fiber_count << " fibers using Construct_Fibers" << std::endl;
528  }
529  else if (get_geom_v3()->get_config() == SpacalGeom_t::kFullProjective_2DTaper_SameLengthFiberPerTower)
530  {
531  int fiber_count = Construct_Fibers_SameLengthFiberPerTower(g_tower,
532  block_logic);
533 
534  if (get_geom_v3()->get_construction_verbose() >= 2)
535  std::cout << "PHG4FullProjSpacalDetector::Construct_Tower::" << GetName()
536  << " - constructed tower ID " << g_tower.id << " with "
537  << fiber_count
538  << " fibers using Construct_Fibers_SameLengthFiberPerTower" << std::endl;
539  }
540  else
541  {
542  G4ExceptionDescription message;
543  message << "can not recognize configuration type " << get_geom_v3()->get_config();
544 
545  G4Exception("PHG4FullProjSpacalDetector::Construct_Tower", "Wrong",
546  FatalException, message, "");
547  }
548 
549  return block_logic;
550 }
551 
552 void PHG4FullProjSpacalDetector::Print(const std::string& /*what*/) const
553 {
554  std::cout << "PHG4FullProjSpacalDetector::Print::" << GetName()
555  << " - Print Geometry:" << std::endl;
556  get_geom_v3()->Print();
557 
558  return;
559 }