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