Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4SpacalPrototype4Detector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4SpacalPrototype4Detector.cc
1 #include "PHG4SpacalPrototype4Detector.h"
2 
3 #include <g4detectors/PHG4CylinderGeomContainer.h>
4 #include <g4detectors/PHG4CylinderGeom_Spacalv1.h> // for PHG4CylinderGeom_Spacalv1:...
5 #include <g4detectors/PHG4CylinderGeom_Spacalv3.h> // for PHG4CylinderGeom_...
6 
7 #include <phparameter/PHParameters.h>
8 
9 #include <g4main/PHG4Detector.h> // for PHG4Detector
10 #include <g4main/PHG4Utils.h>
11 
12 #include <phool/PHCompositeNode.h>
13 #include <phool/PHIODataNode.h>
14 #include <phool/PHNode.h> // for PHNode
15 #include <phool/PHNodeIterator.h> // for PHNodeIterator
16 #include <phool/PHObject.h> // for PHObject
17 #include <phool/getClass.h>
18 
19 #include <Geant4/G4Box.hh>
20 #include <Geant4/G4Colour.hh> // for G4Colour
21 #include <Geant4/G4DisplacedSolid.hh>
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/G4SubtractionSolid.hh>
28 #include <Geant4/G4SystemOfUnits.hh>
29 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
30 #include <Geant4/G4Transform3D.hh> // for G4Transform3D, G4Translate3D
31 #include <Geant4/G4Trap.hh>
32 #include <Geant4/G4Tubs.hh>
33 #include <Geant4/G4Types.hh> // for G4double
34 #include <Geant4/G4UserLimits.hh>
35 #include <Geant4/G4Vector3D.hh>
36 #include <Geant4/G4VisAttributes.hh>
37 
38 #include <boost/foreach.hpp>
39 
40 #include <algorithm>
41 #include <cassert>
42 #include <cmath>
43 #include <cstdlib> // for exit
44 #include <iostream> // for operator<<, basic_ostream
45 #include <limits> // for numeric_limits
46 #include <memory> // for allocator_traits<...
47 #include <numeric> // std::accumulate
48 #include <sstream>
49 #include <string> // std::string, std::to_string
50 #include <vector> // for vector
51 
52 class PHG4CylinderGeom;
53 
54 using namespace std;
55 
56 //_______________________________________________________________
57 //note this inactive thickness is ~1.5% of a radiation length
58 PHG4SpacalPrototype4Detector::PHG4SpacalPrototype4Detector(PHG4Subsystem* subsys, PHCompositeNode* Node, PHParameters* parameters, const std::string& dnam)
59  : PHG4Detector(subsys, Node, dnam)
60  , construction_params(parameters)
61  , cylinder_solid(nullptr)
62  , cylinder_logic(nullptr)
63  , cylinder_physi(nullptr)
64  , //
65  active(0)
66  , absorberactive(0)
67  , //
68  step_limits(nullptr)
69  , clading_step_limits(nullptr)
70  , fiber_core_step_limits(nullptr)
71  , //
72  _geom(nullptr)
73 {
74  SetActive(construction_params->get_int_param("active"));
75  SetAbsorberActive(construction_params->get_int_param("absorberactive"));
76  SuperDetector(dnam);
77 }
78 
79 PHG4SpacalPrototype4Detector::~PHG4SpacalPrototype4Detector(void)
80 {
81  // deleting nullptr pointers is allowed (results in NOOP)
82  // so checking for not null before deleting is not needed
83  delete step_limits;
84  delete clading_step_limits;
85  delete fiber_core_step_limits;
86  delete _geom;
87 }
88 
89 //_______________________________________________________________
90 int PHG4SpacalPrototype4Detector::IsInCylinderActive(
91  const G4VPhysicalVolume* volume)
92 {
93  // cout << "checking detector" << endl;
94  if (active && fiber_core_vol.find(volume) != fiber_core_vol.end())
95  {
96  // return fiber_core_vol.find(volume)->second;
97  return FIBER_CORE;
98  }
99  if (absorberactive)
100  {
101  if (fiber_vol.find(volume) != fiber_vol.end())
102  return FIBER_CLADING;
103 
104  if (block_vol.find(volume) != block_vol.end())
105  return ABSORBER;
106 
107  if (calo_vol.find(volume) != calo_vol.end())
108  return SUPPORT;
109  }
110  return INACTIVE;
111 }
112 
113 //_______________________________________________________________
114 void PHG4SpacalPrototype4Detector::ConstructMe(G4LogicalVolume* logicWorld)
115 {
116  PHNodeIterator iter(topNode());
117  PHCompositeNode* parNode = dynamic_cast<PHCompositeNode*>(iter.findFirst(
118  "PHCompositeNode", "RUN"));
119  assert(parNode);
120 
121  if (!_geom)
122  {
123  _geom = new SpacalGeom_t();
124  }
125  assert(_geom);
126 
127  _geom->set_config(
128  SpacalGeom_t::kFullProjective_2DTaper_Tilted_SameLengthFiberPerTower);
129  // _geom->load_demo_sector_tower_map4();
130  // _geom->set_construction_verbose(2);
131  _geom->ImportParameters(*construction_params);
132 
133  _geom->subtower_consistency_check();
134 
135  // step_limits = new G4UserLimits(_geom->get_calo_step_size() * cm);
136  //
137  // clading_step_limits = new G4UserLimits(
138  // _geom->get_fiber_clading_step_size() * cm);
139  step_limits = nullptr;
140  clading_step_limits = nullptr;
141 
142  fiber_core_step_limits = new G4UserLimits(
143  _geom->get_fiber_core_step_size() * cm);
144 
145  const G4double z_rotation = construction_params->get_double_param(
146  "z_rotation_degree") *
147  degree;
148  const G4double enclosure_x = construction_params->get_double_param(
149  "enclosure_x") *
150  cm;
151  const G4double enclosure_x_shift = construction_params->get_double_param(
152  "enclosure_x_shift") *
153  cm;
154  const G4double enclosure_y = construction_params->get_double_param(
155  "enclosure_y") *
156  cm;
157  const G4double enclosure_z = construction_params->get_double_param(
158  "enclosure_z") *
159  cm;
160 
161  const G4double box_x_shift = (_geom->get_radius() + 0.5 * _geom->get_thickness()) * cm + enclosure_x_shift;
162  const G4double box_z_shift = 0.5 * (_geom->get_zmin() + _geom->get_zmax()) * cm;
163 
164  // G4Tubs* _cylinder_solid = new G4Tubs(G4String(GetName()),
165  // _geom->get_radius() * cm, _geom->get_max_radius() * cm,
166  // _geom->get_length() * cm / 2.0, 0, twopi);
167  // G4VSolid* cylinder_solid = new G4Box(G4String(GetName()),
168  // _geom->get_thickness() * 1.1 * 0.5 * cm, //
169  // twopi / _geom->get_azimuthal_n_sec() * _geom->get_sector_map().size()
170  // * 0.5 * _geom->get_max_radius() * cm, //
171  // _geom->get_length() * cm / 2.0);
172  G4VSolid* cylinder_solid = new G4Box(G4String(GetName()),
173  enclosure_x * 0.5, //
174  enclosure_y * 0.5, //
175  enclosure_z * 0.5);
176 
177  cylinder_solid = new G4DisplacedSolid(G4String(GetName() + "_displaced"),
178  cylinder_solid, 0, //
179  G4ThreeVector(box_x_shift, 0, box_z_shift));
180 
181  G4Material* cylinder_mat = G4Material::GetMaterial("G4_AIR");
182  assert(cylinder_mat);
183 
184  cylinder_logic = new G4LogicalVolume(cylinder_solid, cylinder_mat,
185  G4String(GetName()), 0, 0, 0);
186  G4VisAttributes* VisAtt = new G4VisAttributes();
187  VisAtt->SetColour(G4Colour::White());
188  VisAtt->SetVisibility(true);
189  VisAtt->SetForceSolid(false);
190  VisAtt->SetForceWireframe(true);
191  cylinder_logic->SetVisAttributes(VisAtt);
192 
193  G4Transform3D cylinder_place(
194  G4Translate3D(_geom->get_xpos() * cm, _geom->get_ypos() * cm,
195  _geom->get_zpos() * cm) //
196  * G4RotateZ3D(z_rotation) //
197  * G4Translate3D(
198  -(_geom->get_radius() + 0.5 * _geom->get_thickness()) * cm, 0,
199  0));
200 
201  cylinder_physi = new G4PVPlacement(cylinder_place, cylinder_logic,
202  G4String(GetName()), logicWorld, false, 0, OverlapCheck());
203 
204  // install sectors
205  std::pair<G4LogicalVolume*, G4Transform3D> psec = Construct_AzimuthalSeg();
206  G4LogicalVolume* sec_logic = psec.first;
207  const G4Transform3D& sec_trans = psec.second;
208  BOOST_FOREACH (const SpacalGeom_t::sector_map_t::value_type& val, _geom->get_sector_map())
209  {
210  const int sec = val.first;
211  const double rot = val.second;
212 
213  G4Transform3D sec_place = G4RotateZ3D(rot) * sec_trans;
214 
215  ostringstream name;
216  name << GetName() << "_sec" << sec;
217 
218  G4PVPlacement* calo_phys = new G4PVPlacement(sec_place, sec_logic,
219  G4String(name.str()), cylinder_logic, false, sec,
220  OverlapCheck());
221  calo_vol[calo_phys] = sec;
222  }
223  _geom->set_nscint(_geom->get_nscint() * _geom->get_sector_map().size());
224 
225  // install electronics
226  const G4double electronics_thickness = construction_params->get_double_param(
227  "electronics_thickness") *
228  cm;
229  const string electronics_material = construction_params->get_string_param(
230  "electronics_material");
231 
232  G4VSolid* electronics_solid = new G4Box(G4String(GetName()),
233  electronics_thickness / 2.0, //
234  sin(
235  twopi / _geom->get_azimuthal_n_sec() * _geom->get_sector_map().size() / 2) *
236  (_geom->get_radius() - _geom->get_max_lightguide_height()) * cm, //
237  _geom->get_length() * cm / 2.0);
238 
239  electronics_solid = new G4DisplacedSolid(G4String(GetName() + "_displaced"),
240  electronics_solid,
241  0, //
242  G4ThreeVector(
243  cos(
244  twopi / _geom->get_azimuthal_n_sec() * _geom->get_sector_map().size() / 2) *
245  (_geom->get_radius() - _geom->get_max_lightguide_height()) * cm -
246  electronics_thickness,
247  0, box_z_shift));
248 
249  G4Material* electronics_mat = G4Material::GetMaterial(electronics_material);
250  assert(electronics_mat);
251 
252  G4LogicalVolume* electronics_logic = new G4LogicalVolume(electronics_solid,
253  electronics_mat, G4String(GetName()) + "_Electronics", 0, 0, 0);
254  VisAtt = new G4VisAttributes();
255  PHG4Utils::SetColour(VisAtt, electronics_material);
256  VisAtt->SetVisibility(true);
257  VisAtt->SetForceSolid(not _geom->is_virualize_fiber());
258  electronics_logic->SetVisAttributes(VisAtt);
259 
260  new G4PVPlacement(G4Transform3D::Identity, electronics_logic,
261  G4String(GetName()) + "_Electronics", cylinder_logic, false, 0,
262  OverlapCheck());
263 
264  // install the enclosure box
265  const G4double enclosure_thickness = construction_params->get_double_param(
266  "enclosure_thickness") *
267  cm;
268  const string enclosure_material = construction_params->get_string_param(
269  "enclosure_material");
270 
271  G4VSolid* outer_enclosur_solid = new G4Box(
272  G4String(GetName()) + "_outer_enclosur_solid", enclosure_x * 0.5, //
273  enclosure_y * 0.5, //
274  enclosure_z * 0.5);
275  G4VSolid* inner_enclosur_solid = new G4Box(
276  G4String(GetName()) + "_inner_enclosur_solid",
277  enclosure_x * 0.5 - enclosure_thickness, //
278  enclosure_y * 0.5 - enclosure_thickness, //
279  enclosure_z * 0.5 - enclosure_thickness);
280  G4VSolid* enclosure_solid = new G4SubtractionSolid(
281  G4String(GetName()) + "_enclosure_solid", //
282  outer_enclosur_solid, inner_enclosur_solid);
283  enclosure_solid = new G4DisplacedSolid(
284  G4String(GetName() + "_enclosure_solid_displaced"), enclosure_solid, 0, //
285  G4ThreeVector(box_x_shift, 0, box_z_shift));
286 
287  G4Material* enclosure_mat = G4Material::GetMaterial(enclosure_material);
288  assert(enclosure_mat);
289 
290  G4LogicalVolume* enclosure_logic = new G4LogicalVolume(enclosure_solid,
291  enclosure_mat, G4String(GetName()) + "_enclosure", 0, 0, 0);
292  VisAtt = new G4VisAttributes();
293  PHG4Utils::SetColour(VisAtt, enclosure_material);
294  VisAtt->SetVisibility(true);
295  VisAtt->SetForceSolid(not _geom->is_virualize_fiber());
296  enclosure_logic->SetVisAttributes(VisAtt);
297 
298  new G4PVPlacement(G4Transform3D::Identity, enclosure_logic,
299  G4String(GetName()) + "_enclosure", cylinder_logic, false, 0,
300  OverlapCheck());
301 
302  if (active)
303  {
304  ostringstream geonode;
305  if (superdetector != "NONE")
306  {
307  geonode << "CYLINDERGEOM_" << superdetector;
308  }
309  else
310  {
311  geonode << "CYLINDERGEOM_" << detector_type;
312  }
313  PHG4CylinderGeomContainer* geo = findNode::getClass<
314  PHG4CylinderGeomContainer>(topNode(), geonode.str());
315  if (!geo)
316  {
317  geo = new PHG4CylinderGeomContainer();
318  PHNodeIterator iter(topNode());
319  PHCompositeNode* runNode =
320  dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode",
321  "RUN"));
322  PHIODataNode<PHObject>* newNode = new PHIODataNode<PHObject>(geo,
323  geonode.str(), "PHObject");
324  runNode->addNode(newNode);
325  }
326  // here in the detector class we have internal units, convert to cm
327  // before putting into the geom object
328  PHG4CylinderGeom* mygeom = new SpacalGeom_t(*_geom);
329  geo->AddLayerGeom(0, mygeom);
330  if (_geom->get_construction_verbose() >= 1)
331  {
332  cout << "PHG4SpacalPrototype4Detector::Construct::" << GetName()
333  << " - Print Layer Geometry:" << endl;
334  geo->identify();
335  }
336  }
337 
338  if (absorberactive)
339  {
340  ostringstream geonode;
341  if (superdetector != "NONE")
342  {
343  geonode << "CYLINDERGEOM_ABSORBER_" << superdetector;
344  }
345  else
346  {
347  geonode << "CYLINDERGEOM_ABSORBER_" << detector_type << "_" << 0;
348  }
349  PHG4CylinderGeomContainer* geo = findNode::getClass<PHG4CylinderGeomContainer>(topNode(), geonode.str());
350  if (!geo)
351  {
352  geo = new PHG4CylinderGeomContainer();
353  PHNodeIterator iter(topNode());
354  PHCompositeNode* runNode =
355  dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode",
356  "RUN"));
357  PHIODataNode<PHObject>* newNode = new PHIODataNode<PHObject>(geo, geonode.str(), "PHObject");
358  runNode->addNode(newNode);
359  }
360  // here in the detector class we have internal units, convert to cm
361  // before putting into the geom object
362  PHG4CylinderGeom* mygeom = new SpacalGeom_t(*_geom);
363  geo->AddLayerGeom(0, mygeom);
364  if (_geom->get_construction_verbose() >= 1)
365  {
366  cout << "PHG4SpacalPrototype4Detector::Construct::" << GetName()
367  << " - Print Absorber Layer Geometry:" << endl;
368  geo->identify();
369  }
370  }
371 
372  if (_geom->get_construction_verbose() >= 1)
373  {
374  cout << "PHG4SpacalPrototype4Detector::Construct::" << GetName()
375  << " - Completed. Print G4 Geometry:" << endl;
376  Print();
377  cout << "box_x_shift = " << box_x_shift << endl;
378  cout << "box_z_shift = " << box_z_shift << endl;
379  }
380 }
381 
382 std::pair<G4LogicalVolume*, G4Transform3D>
383 PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg()
384 {
385  assert(_geom);
386  assert(_geom->get_azimuthal_n_sec() > 4);
387 
388  // basic tilt geometry
389  const G4double half_chord_backend =
390  _geom->get_max_radius() * cm * tan(pi / _geom->get_azimuthal_n_sec()) //
391  + fabs(_geom->get_thickness() * cm * 0.5 * tan(_geom->get_azimuthal_tilt()));
392  const G4double reduced_outer_radius = sqrt(pow(_geom->get_max_radius() * cm, 2) - half_chord_backend * half_chord_backend);
393  const G4double enclosure_depth = reduced_outer_radius - _geom->get_radius() * cm;
394  const G4double enclosure_center = 0.5 * (reduced_outer_radius + _geom->get_radius() * cm);
395  const G4double enclosure_half_height_half_width = enclosure_center * tan(pi / _geom->get_azimuthal_n_sec());
396 
397  const G4double width_adj1 = tan(_geom->get_azimuthal_tilt() - pi / _geom->get_azimuthal_n_sec()) * enclosure_depth * 0.5;
398  const G4double width_adj2 = tan(_geom->get_azimuthal_tilt() + pi / _geom->get_azimuthal_n_sec()) * enclosure_depth * 0.5;
399 
400  const G4double center_adj = (width_adj1 + width_adj2) * 0.5;
401  const G4double center_tilt_angle = atan2(center_adj, enclosure_depth * 0.5);
402  const G4double inner_half_width = enclosure_half_height_half_width + 0.5 * (width_adj1 - width_adj2);
403  const G4double outter_half_width = enclosure_half_height_half_width + 0.5 * (-width_adj1 + width_adj2);
404 
405  // enclosure walls
406  const G4double edge1_tilt_angle = atan2(width_adj1, enclosure_depth * 0.5);
407  const G4double edge2_tilt_angle = atan2(width_adj2, enclosure_depth * 0.5);
408  // const G4double edge1_half_depth = sqrt(width_adj1 * width_adj1 + enclosure_depth * enclosure_depth * .25);
409  // const G4double edge2_half_depth = sqrt(width_adj2 * width_adj2 + enclosure_depth * enclosure_depth * .25);
410 
411  // projective center
412  const G4double half_projection_ratio = 0.5 * (-width_adj1 + width_adj2) / enclosure_half_height_half_width;
413  const G4double projection_center_y = enclosure_center - ((enclosure_depth * 0.5) / half_projection_ratio);
414  const G4double projection_center_x = center_adj / half_projection_ratio;
415 
416  // blocks azimuthal segmentation
417  const int phi_bin_in_sec = _geom->get_max_phi_bin_in_sec();
418  assert(phi_bin_in_sec >= 1);
419  const G4double block_azimuth_angle = (edge2_tilt_angle - edge1_tilt_angle) / phi_bin_in_sec;
420  assert(block_azimuth_angle > 0);
421  if (!(fabs(block_azimuth_angle - M_PI * 2 / _geom->get_azimuthal_n_sec() / phi_bin_in_sec) < M_PI * numeric_limits<G4double>::epsilon()))
422  {
423  cout << "angle/nsec out of range: " << M_PI * numeric_limits<G4double>::epsilon() << endl;
424  exit(1);
425  }
426  const G4double block_edge1_half_width = enclosure_half_height_half_width - (_geom->get_sidewall_thickness() * cm + _geom->get_sidewall_outer_torr() * cm + 2.0 * _geom->get_assembly_spacing() * cm) / cos(edge1_tilt_angle);
427  const G4double block_edge2_half_width = enclosure_half_height_half_width - (_geom->get_sidewall_thickness() * cm + _geom->get_sidewall_outer_torr() * cm + 2.0 * _geom->get_assembly_spacing() * cm) / cos(edge2_tilt_angle);
428  G4double block_width_ratio = 0;
429  for (int s = 0; s < phi_bin_in_sec; ++s)
430  {
431  block_width_ratio += 1 / cos(block_azimuth_angle * (0.5 + s) + edge1_tilt_angle);
432  }
433  const G4double block_half_height_width = (block_edge1_half_width + block_edge2_half_width) / block_width_ratio;
434  assert(block_half_height_width > 0);
435 
436  // write out the azimuthal block geometry
437  // block azimuth geometry records
438  struct block_azimuth_geom
439  {
440  G4double angle;
441  G4double projection_center_y;
442  G4double projection_center_x;
443  G4double projection_length;
444  };
445  vector<block_azimuth_geom> block_azimuth_geoms(phi_bin_in_sec,
446  block_azimuth_geom{
447  numeric_limits<double>::signaling_NaN(),
448  numeric_limits<double>::signaling_NaN(),
449  numeric_limits<double>::signaling_NaN(),
450  numeric_limits<double>::signaling_NaN()}); // [phi-bin in sector] -> azimuth geometry
451  G4double block_x_edge1 = block_edge1_half_width;
452  for (int s = 0; s < phi_bin_in_sec; ++s)
453  {
454  block_azimuth_geom& geom = block_azimuth_geoms[s];
455 
456  geom.angle = block_azimuth_angle * (0.5 + s) + edge1_tilt_angle;
457  const G4double block_x_size = block_half_height_width / cos(geom.angle);
458  assert(block_x_size > 0);
459  const G4double x_center = block_x_edge1 - 0.5 * block_x_size;
460 
461  // projection center per block
462  geom.projection_length = block_half_height_width / 2. / tan(block_azimuth_angle / 2.);
463  assert(geom.projection_length > 0);
464  geom.projection_center_y = enclosure_center - geom.projection_length * cos(geom.angle);
465  geom.projection_center_x = x_center + geom.projection_length * sin(geom.angle);
466 
467  // next step
468  block_x_edge1 -= block_x_size;
469  }
470 
471  //write out the azimuthal block divider's geometry
472  struct block_divider_azimuth_geom
473  {
474  G4double angle;
475  G4double projection_center_y;
476  G4double projection_center_x;
477  G4double thickness; // thickness in the approximate azimuth direction
478  G4double radial_displacement;
479  G4double width;
480  };
481  assert(phi_bin_in_sec >= 1);
482  vector<block_divider_azimuth_geom> divider_azimuth_geoms(phi_bin_in_sec - 1,
483  block_divider_azimuth_geom{
484  numeric_limits<double>::signaling_NaN(),
485  numeric_limits<double>::signaling_NaN(),
486  numeric_limits<double>::signaling_NaN(),
487  numeric_limits<double>::signaling_NaN(),
488  numeric_limits<double>::signaling_NaN(),
489  numeric_limits<double>::signaling_NaN()});
490 
491  if (_geom->get_sidewall_thickness() > 0)
492  {
493  for (int s = 0; s < phi_bin_in_sec - 1; ++s)
494  {
495  block_divider_azimuth_geom& geom = divider_azimuth_geoms[s];
496 
497  geom.angle = 0.5 * (block_azimuth_geoms[s].angle + block_azimuth_geoms[s + 1].angle);
498  geom.projection_center_y = 0.5 * (block_azimuth_geoms[s].projection_center_y + block_azimuth_geoms[s + 1].projection_center_y);
499  geom.projection_center_x = 0.5 * (block_azimuth_geoms[s].projection_center_x + block_azimuth_geoms[s + 1].projection_center_x);
500  geom.radial_displacement = 0.5 * (block_azimuth_geoms[s].projection_length + block_azimuth_geoms[s + 1].projection_length);
501 
502  geom.thickness = 2.0 * _geom->get_assembly_spacing() * cm * cos(block_azimuth_angle / 2.) - 2 * um;
503  geom.width = _geom->get_divider_width() * cm;
504  }
505  }
506 
507  if (fabs(block_x_edge1 - (-block_edge2_half_width)) > _geom->get_assembly_spacing() * cm)
508  {
509  cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg - ERROR - " << endl
510  << "\t block_x_edge1 = " << block_x_edge1 << endl
511  << "\t block_edge2_half_width = " << block_edge2_half_width << endl
512  << "\t fabs(block_x_edge1 - (-block_edge2_half_width)) = " << fabs(block_x_edge1 - (-block_edge2_half_width)) << endl
513  << "\t _geom->get_assembly_spacing() * cm = " << _geom->get_assembly_spacing() * cm << endl;
514  }
515  if (!(fabs(block_x_edge1 - (-block_edge2_half_width)) < _geom->get_assembly_spacing() * cm)) // closure check
516  {
517  cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg - ERROR - closure check failed: " << fabs(block_x_edge1 - (-block_edge2_half_width)) << endl;
518  exit(1);
519  }
520 
521  if (Verbosity())
522  {
523  cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg - " << endl
524  << "\t edge1_tilt_angle = " << edge1_tilt_angle << endl
525  << "\t edge2_tilt_angle = " << edge2_tilt_angle << endl
526  << "\t projection_center_y = " << projection_center_y << endl
527  << "\t projection_center_x = " << projection_center_x << endl
528  << "\t block_azimuth_angle = " << block_azimuth_angle << endl
529  << "\t block_edge1_half_width = " << block_edge1_half_width << endl
530  << "\t block_edge2_half_width = " << block_edge2_half_width << endl
531  << "\t block_width_ratio = " << block_width_ratio << endl
532  << "\t block_half_height_width = " << block_half_height_width << endl;
533 
534  for (int s = 0; s < phi_bin_in_sec; ++s)
535  {
536  cout << "\t block[" << s << "].angle = " << block_azimuth_geoms[s].angle << endl;
537  cout << "\t block[" << s << "].projection_center_y = " << block_azimuth_geoms[s].projection_center_y << endl;
538  cout << "\t block[" << s << "].projection_center_x = " << block_azimuth_geoms[s].projection_center_x << endl;
539  }
540  for (int s = 0; s < phi_bin_in_sec - 1; ++s)
541  {
542  cout << "\t divider[" << s << "].angle = " << divider_azimuth_geoms[s].angle << endl;
543  cout << "\t divider[" << s << "].projection_center_x = " << divider_azimuth_geoms[s].projection_center_x << endl;
544  cout << "\t divider[" << s << "].projection_center_y = " << divider_azimuth_geoms[s].projection_center_y << endl;
545  cout << "\t divider[" << s << "].radial_displacement = " << divider_azimuth_geoms[s].radial_displacement << endl;
546  cout << "\t divider[" << s << "].thickness = " << divider_azimuth_geoms[s].thickness << endl;
547  cout << "\t divider[" << s << "].width = " << divider_azimuth_geoms[s].width << endl;
548  }
549  }
550 
551  assert(enclosure_depth > 10 * cm);
552 
553  G4VSolid* sec_solid = new G4Trap(
554  /*const G4String& pName*/ G4String(GetName() + string("_sec_trap")),
555  enclosure_depth * 0.5, // G4double pDz,
556  center_tilt_angle, halfpi, // G4double pTheta, G4double pPhi,
557  inner_half_width, _geom->get_length() * cm / 2.0, _geom->get_length() * cm / 2.0, // G4double pDy1, G4double pDx1, G4double pDx2,
558  0, // G4double pAlp1,
559  outter_half_width, _geom->get_length() * cm / 2.0, _geom->get_length() * cm / 2.0, // G4double pDy2, G4double pDx3, G4double pDx4,
560  0 // G4double pAlp2 //
561  );
562  const G4double box_z_shift = 0.5 * (_geom->get_zmin() + _geom->get_zmax()) * cm;
563  G4Transform3D sec_solid_transform =
564  G4TranslateZ3D(box_z_shift) * G4TranslateY3D(enclosure_center) * G4RotateY3D(halfpi) * G4RotateX3D(-halfpi);
565  G4VSolid* sec_solid_place = new G4DisplacedSolid(
566  G4String(GetName() + string("_sec")), sec_solid, sec_solid_transform);
567 
568  G4Material* cylinder_mat = G4Material::GetMaterial("G4_AIR");
569  assert(cylinder_mat);
570 
571  G4LogicalVolume* sec_logic = new G4LogicalVolume(sec_solid_place, cylinder_mat,
572  G4String(G4String(GetName() + string("_sec"))), 0, 0, nullptr);
573 
574  G4VisAttributes* VisAtt = new G4VisAttributes();
575  VisAtt->SetColor(.5, .9, .5, .1);
576  VisAtt->SetVisibility(true);
577  VisAtt->SetForceSolid(false);
578  VisAtt->SetForceWireframe(true);
579  sec_logic->SetVisAttributes(VisAtt);
580 
581  // // construct walls
582  //
583  // G4Material* wall_mat = G4Material::GetMaterial(_geom->get_sidewall_mat());
584  // assert(wall_mat);
585  //
586  // if (_geom->get_sidewall_thickness() > 0)
587  // {
588  // // end walls
589  // if (_geom->get_construction_verbose() >= 1)
590  // {
591  // cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg::" << GetName()
592  // << " - construct end walls." << endl;
593  // }
594  // // G4Tubs* wall_solid = new G4Tubs(G4String(GetName() + string("_EndWall")),
595  // // _geom->get_radius() * cm + _geom->get_sidewall_outer_torr() * cm,
596  // // _geom->get_max_radius() * cm - _geom->get_sidewall_outer_torr() * cm,
597  // // _geom->get_sidewall_thickness() * cm / 2.0,
598  // // halfpi - pi / _geom->get_azimuthal_n_sec(),
599  // // twopi / _geom->get_azimuthal_n_sec());
600  // const G4double side_wall_half_thickness = _geom->get_sidewall_thickness() * cm / 2.0;
601  // G4VSolid* wall_solid = new G4Trap(G4String(GetName() + string("_EndWall_trap")),
602  // enclosure_depth * 0.5, // G4double pDz,
603  // center_tilt_angle, halfpi, // G4double pTheta, G4double pPhi,
604  // inner_half_width, side_wall_half_thickness, side_wall_half_thickness, // G4double pDy1, G4double pDx1, G4double pDx2,
605  // 0, // G4double pAlp1,
606  // outter_half_width, side_wall_half_thickness, side_wall_half_thickness, // G4double pDy2, G4double pDx3, G4double pDx4,
607  // 0 // G4double pAlp2 //
608  // );
609  // G4VSolid* wall_solid_place = new G4DisplacedSolid(
610  // G4String(GetName() + string("_EndWall")), wall_solid, sec_solid_transform);
611  //
612  // G4LogicalVolume* wall_logic = new G4LogicalVolume(wall_solid_place, wall_mat,
613  // G4String(G4String(GetName() + string("_EndWall"))), 0, 0,
614  // nullptr);
615  // GetDisplayAction()->AddVolume(wall_logic, "Wall");
616  //
617  // typedef map<int, double> z_locations_t;
618  // z_locations_t z_locations;
619  // z_locations[000] = _geom->get_sidewall_thickness() * cm / 2.0 + _geom->get_assembly_spacing() * cm;
620  // z_locations[001] = _geom->get_length() * cm / 2.0 - (_geom->get_sidewall_thickness() * cm / 2.0 + _geom->get_assembly_spacing() * cm);
621  // z_locations[100] = -(_geom->get_sidewall_thickness() * cm / 2.0 + _geom->get_assembly_spacing() * cm);
622  // z_locations[101] = -(_geom->get_length() * cm / 2.0 - (_geom->get_sidewall_thickness() * cm / 2.0 + _geom->get_assembly_spacing() * cm));
623  //
624  // BOOST_FOREACH (z_locations_t::value_type& val, z_locations)
625  // {
626  // if (_geom->get_construction_verbose() >= 2)
627  // cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg::"
628  // << GetName() << " - constructed End Wall ID " << val.first
629  // << " @ Z = " << val.second << endl;
630  //
631  // G4Transform3D wall_trans = G4TranslateZ3D(val.second);
632  //
633  // G4PVPlacement* wall_phys = new G4PVPlacement(wall_trans, wall_logic,
634  // G4String(GetName()) + G4String("_EndWall_") + to_string(val.first), sec_logic,
635  // false, val.first, OverlapCheck());
636  //
637  // calo_vol[wall_phys] = val.first;
638  // assert(gdml_config);
639  // gdml_config->exclude_physical_vol(wall_phys);
640  // }
641  // }
642  // //
643  // if (_geom->get_sidewall_thickness() > 0)
644  // {
645  // // side walls
646  // if (_geom->get_construction_verbose() >= 1)
647  // {
648  // cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg::" << GetName()
649  // << " - construct side walls." << endl;
650  // }
651  //
652  // typedef map<int, pair<int, int> > sign_t;
653  // sign_t signs;
654  // signs[100] = make_pair(+1, +1);
655  // signs[101] = make_pair(+1, -1);
656  // signs[200] = make_pair(-1, +1);
657  // signs[201] = make_pair(-1, -1);
658  //
659  // BOOST_FOREACH (sign_t::value_type& val, signs)
660  // {
661  // const int sign_z = val.second.first;
662  // const int sign_azimuth = val.second.second;
663  //
664  // if (_geom->get_construction_verbose() >= 2)
665  // cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg::"
666  // << GetName() << " - constructed Side Wall ID " << val.first
667  // << " with"
668  // << " Shift X = "
669  // << sign_azimuth * (_geom->get_sidewall_thickness() * cm / 2.0 + _geom->get_sidewall_outer_torr() * cm)
670  // << " Rotation Z = "
671  // << sign_azimuth * pi / _geom->get_azimuthal_n_sec()
672  // << " Shift Z = " << sign_z * (_geom->get_length() * cm / 4)
673  // << endl;
674  //
675  // const G4double azimuth_roate =
676  // sign_azimuth > 0 ? edge1_tilt_angle : edge2_tilt_angle;
677  // const G4double edge_half_depth = -_geom->get_sidewall_thickness() * cm - _geom->get_sidewall_outer_torr() * cm +
678  // (sign_azimuth > 0 ? edge1_half_depth : edge2_half_depth);
679  //
680  // G4Box* wall_solid = new G4Box(G4String(GetName() + G4String("_SideWall_") + to_string(val.first)),
681  // _geom->get_sidewall_thickness() * cm / 2.0,
682  // edge_half_depth,
683  // (_geom->get_length() / 2. - 2 * (_geom->get_sidewall_thickness() + 2. * _geom->get_assembly_spacing())) * cm * .5);
684  //
685  // G4LogicalVolume* wall_logic = new G4LogicalVolume(wall_solid, wall_mat,
686  // G4String(G4String(GetName() + G4String("_SideWall_") + to_string(val.first))), 0, 0,
687  // nullptr);
688  // GetDisplayAction()->AddVolume(wall_logic, "Wall");
689  //
690  // const G4Transform3D wall_trans =
691  // G4TranslateZ3D(sign_z * (_geom->get_length() * cm / 4)) *
692  // G4TranslateY3D(enclosure_center) *
693  // G4TranslateX3D(sign_azimuth * enclosure_half_height_half_width) *
694  // G4RotateZ3D(azimuth_roate) *
695  // G4TranslateX3D(-sign_azimuth * (_geom->get_sidewall_thickness() * cm / 2.0 + _geom->get_sidewall_outer_torr() * cm));
696  //
697  // G4PVPlacement* wall_phys = new G4PVPlacement(wall_trans, wall_logic,
698  // G4String(GetName()) + G4String("_SideWall_") + to_string(val.first), sec_logic,
699  // false, val.first, OverlapCheck());
700  //
701  // calo_vol[wall_phys] = val.first;
702  // assert(gdml_config);
703  // gdml_config->exclude_physical_vol(wall_phys);
704  // }
705  // }
706  //
707  // // construct dividers
708  // if (_geom->get_divider_width() > 0)
709  // {
710  // if (_geom->get_construction_verbose() >= 1)
711  // {
712  // cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg::" << GetName()
713  // << " - construct dividers" << endl;
714  // }
715  //
716  // G4Material* divider_mat = G4Material::GetMaterial(_geom->get_divider_mat());
717  // assert(divider_mat);
718  //
719  // int ID = 300;
720  // for (const auto& geom : divider_azimuth_geoms)
721  // {
722  // G4Box* divider_solid = new G4Box(G4String(GetName() + G4String("_Divider_") + to_string(ID)),
723  // geom.thickness / 2.0,
724  // geom.width / 2.,
725  // (_geom->get_length() / 2. - 2 * (_geom->get_sidewall_thickness() + 2. * _geom->get_assembly_spacing())) * cm * .5);
726  //
727  // G4LogicalVolume* wall_logic = new G4LogicalVolume(divider_solid, divider_mat,
728  // G4String(G4String(GetName() + G4String("_Divider_") + to_string(ID))), 0, 0,
729  // nullptr);
730  // GetDisplayAction()->AddVolume(wall_logic, "Divider");
731  //
732  // for (int sign_z = -1; sign_z <= 1; sign_z += 2)
733  // {
734  // G4Transform3D wall_trans =
735  // G4TranslateX3D(geom.projection_center_x) *
736  // G4TranslateY3D(geom.projection_center_y) *
737  // G4RotateZ3D(geom.angle) *
738  // G4TranslateY3D(geom.radial_displacement) *
739  // G4TranslateZ3D(sign_z * (_geom->get_length() * cm / 4));
740  //
741  // G4PVPlacement* wall_phys = new G4PVPlacement(wall_trans, wall_logic,
742  // G4String(GetName()) + G4String("_Divider_") + to_string(ID), sec_logic,
743  // false, ID, OverlapCheck());
744  //
745  // calo_vol[wall_phys] = ID;
746  // // assert(gdml_config);
747  // // gdml_config->exclude_physical_vol(wall_phys);
748  //
749  // if (Verbosity())
750  // {
751  // cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg - placing divider " << wall_phys->GetName() << " copy ID " << ID << endl;
752  // }
753  //
754  // ++ID;
755  // }
756  // } // for (const auto & geom : divider_azimuth_geoms)
757  // }
758 
759  // // construct towers
760  //
761  BOOST_FOREACH (const SpacalGeom_t::tower_map_t::value_type& val, _geom->get_sector_tower_map())
762  {
763  SpacalGeom_t::geom_tower g_tower = val.second;
764 
765  const int tower_id = g_tower.id;
766  const int tower_phi_id_in_sec = tower_id % 10;
767  assert(tower_phi_id_in_sec >= 0);
768  assert(tower_phi_id_in_sec < phi_bin_in_sec);
769 
770  const auto& block_azimuth_geom = block_azimuth_geoms.at(tower_phi_id_in_sec);
771 
772  G4LogicalVolume* LV_tower = Construct_Tower(g_tower);
773 
774  G4Transform3D block_trans =
775  G4TranslateX3D(block_azimuth_geom.projection_center_x) *
776  G4TranslateY3D(block_azimuth_geom.projection_center_y) *
777  G4RotateZ3D(block_azimuth_geom.angle) *
778  G4TranslateX3D(g_tower.centralX * cm) *
779  G4TranslateY3D(g_tower.centralY * cm) *
780  G4TranslateZ3D(g_tower.centralZ * cm) *
781  G4RotateX3D(g_tower.pRotationAngleX * rad);
782 
783  const bool overlapcheck_block = OverlapCheck() and (_geom->get_construction_verbose() >= 2);
784 
785  G4PVPlacement* block_phys = new G4PVPlacement(block_trans, LV_tower,
786  G4String(GetName()) + G4String("_Tower_") + to_string(g_tower.id), sec_logic, false,
787  g_tower.id, overlapcheck_block);
788  block_vol[block_phys] = g_tower.id;
789  // assert(gdml_config);
790  // gdml_config->exclude_physical_vol(block_phys);
791 
792  if (g_tower.LightguideHeight > 0)
793  {
794  // also build a light guide
795 
796  for (int ix = 0; ix < g_tower.NSubtowerX; ix++)
797  // int ix = 0;
798  {
799  for (int iy = 0; iy < g_tower.NSubtowerY; iy++)
800  // int iy = 0;
801  {
802  G4LogicalVolume* LV_lg = Construct_LightGuide(g_tower, ix,
803  iy);
804 
805  G4PVPlacement* lg_phys = new G4PVPlacement(block_trans, LV_lg, LV_lg->GetName(),
806  sec_logic, false, g_tower.id, overlapcheck_block);
807 
808  block_vol[lg_phys] = g_tower.id * 100 + ix * 10 + iy;
809 
810  // assert(gdml_config);
811  // gdml_config->exclude_physical_vol(lg_phys);
812  }
813  }
814  }
815  }
816 
817  cout << "PHG4SpacalPrototype4Detector::Construct_AzimuthalSeg::" << GetName()
818  << " - constructed " << _geom->get_sector_tower_map().size()
819  << " unique towers" << endl;
820 
821  return make_pair(sec_logic, G4Transform3D::Identity);
822 }
823 
824 G4LogicalVolume*
825 PHG4SpacalPrototype4Detector::Construct_Fiber(const G4double length,
826  const string& id)
827 {
828  G4Tubs* fiber_solid = new G4Tubs(G4String(GetName() + string("_fiber") + id),
829  0, _geom->get_fiber_outer_r() * cm, length / 2.0, 0, twopi);
830 
831  G4Material* clading_mat = G4Material::GetMaterial(
832  _geom->get_fiber_clading_mat());
833  assert(clading_mat);
834 
835  G4LogicalVolume* fiber_logic = new G4LogicalVolume(fiber_solid, clading_mat,
836  G4String(G4String(GetName() + string("_fiber") + id)), 0, 0,
837  clading_step_limits);
838 
839  {
840  G4VisAttributes* VisAtt = new G4VisAttributes();
841  PHG4Utils::SetColour(VisAtt, "G4_POLYSTYRENE");
842  VisAtt->SetVisibility(_geom->is_virualize_fiber());
843  VisAtt->SetForceSolid(_geom->is_virualize_fiber());
844  fiber_logic->SetVisAttributes(VisAtt);
845  }
846 
847  G4Tubs* core_solid = new G4Tubs(
848  G4String(GetName() + string("_fiber_core") + id), 0,
849  _geom->get_fiber_core_diameter() * cm / 2, length / 2.0, 0, twopi);
850 
851  G4Material* core_mat = G4Material::GetMaterial(_geom->get_fiber_core_mat());
852  assert(core_mat);
853 
854  G4LogicalVolume* core_logic = new G4LogicalVolume(core_solid, core_mat,
855  G4String(G4String(GetName() + string("_fiber_core") + id)), 0, 0,
856  fiber_core_step_limits);
857 
858  {
859  G4VisAttributes* VisAtt = new G4VisAttributes();
860  PHG4Utils::SetColour(VisAtt, "G4_POLYSTYRENE");
861  VisAtt->SetVisibility(false);
862  VisAtt->SetForceSolid(false);
863  core_logic->SetVisAttributes(VisAtt);
864  }
865 
866  const bool overlapcheck_fiber = OverlapCheck() and (_geom->get_construction_verbose() >= 3);
867  G4PVPlacement* core_physi = new G4PVPlacement(0, G4ThreeVector(), core_logic,
868  G4String(G4String(GetName() + string("_fiber_core") + id)), fiber_logic,
869  false, 0, overlapcheck_fiber);
870  fiber_core_vol[core_physi] = 0;
871 
872  return fiber_logic;
873 }
874 
875 void PHG4SpacalPrototype4Detector::Print(const std::string& /*what*/) const
876 {
877  assert(_geom);
878 
879  cout << "PHG4SpacalPrototype4Detector::Print::" << GetName()
880  << " - Print Geometry:" << endl;
881  _geom->Print();
882 
883  return;
884 }
885 
887 int PHG4SpacalPrototype4Detector::Construct_Fibers_SameLengthFiberPerTower(
888  const PHG4SpacalPrototype4Detector::SpacalGeom_t::geom_tower& g_tower,
889  G4LogicalVolume* LV_tower)
890 {
891  assert(_geom);
892 
893  // construct fibers
894 
895  // first check out the fibers geometry
896 
897  typedef map<int, pair<G4Vector3D, G4Vector3D> > fiber_par_map;
898  fiber_par_map fiber_par;
899  G4double min_fiber_length = g_tower.pDz * cm * 4;
900 
901  G4Vector3D v_zshift = G4Vector3D(tan(g_tower.pTheta) * cos(g_tower.pPhi),
902  tan(g_tower.pTheta) * sin(g_tower.pPhi), 1) *
903  g_tower.pDz;
904  // int fiber_ID = 0;
905  for (int ix = 0; ix < g_tower.NFiberX; ix++)
906  // int ix = 0;
907  {
908  const double weighted_ix = static_cast<double>(ix) / (g_tower.NFiberX - 1.);
909 
910  const double weighted_pDx1 = (g_tower.pDx1 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
911  const double weighted_pDx2 = (g_tower.pDx2 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
912 
913  const double weighted_pDx3 = (g_tower.pDx3 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
914  const double weighted_pDx4 = (g_tower.pDx4 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
915 
916  for (int iy = 0; iy < g_tower.NFiberY; iy++)
917  // int iy = 0;
918  {
919  if ((ix + iy) % 2 == 1)
920  continue; // make a triangle pattern
921 
922  const double weighted_iy = static_cast<double>(iy) / (g_tower.NFiberY - 1.);
923 
924  const double weighted_pDy1 = (g_tower.pDy1 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
925  const double weighted_pDy2 = (g_tower.pDy2 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
926 
927  const double weighted_pDx12 = weighted_pDx1 * (1 - weighted_iy) + weighted_pDx2 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp1);
928  const double weighted_pDx34 = weighted_pDx3 * (1 - weighted_iy) + weighted_pDx4 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp2);
929 
930  G4Vector3D v1 = G4Vector3D(weighted_pDx12, weighted_pDy1, 0) - v_zshift;
931  G4Vector3D v2 = G4Vector3D(weighted_pDx34, weighted_pDy2, 0) + v_zshift;
932 
933  G4Vector3D vector_fiber = (v2 - v1);
934  vector_fiber *= (vector_fiber.mag() - _geom->get_fiber_outer_r()) / vector_fiber.mag(); // shrink by fiber boundary protection
935  G4Vector3D center_fiber = (v2 + v1) / 2;
936 
937  // convert to Geant4 units
938  vector_fiber *= cm;
939  center_fiber *= cm;
940 
941  const int fiber_ID = g_tower.compose_fiber_id(ix, iy);
942  fiber_par[fiber_ID] = make_pair(vector_fiber, center_fiber);
943 
944  const G4double fiber_length = vector_fiber.mag();
945 
946  min_fiber_length = min(fiber_length, min_fiber_length);
947 
948  // ++fiber_ID;
949  }
950  }
951 
952  int fiber_count = 0;
953 
954  const G4double fiber_length = min_fiber_length;
955  vector<G4double> fiber_cut;
956 
957  ostringstream ss;
958  ss << string("_Tower") << g_tower.id;
959  G4LogicalVolume* fiber_logic = Construct_Fiber(fiber_length, ss.str());
960 
961  BOOST_FOREACH (const fiber_par_map::value_type& val, fiber_par)
962  {
963  const int fiber_ID = val.first;
964  G4Vector3D vector_fiber = val.second.first;
965  G4Vector3D center_fiber = val.second.second;
966  const G4double optimal_fiber_length = vector_fiber.mag();
967 
968  const G4Vector3D v1 = center_fiber - 0.5 * vector_fiber;
969 
970  // keep a statistics
971  assert(optimal_fiber_length - fiber_length >= 0);
972  fiber_cut.push_back(optimal_fiber_length - fiber_length);
973 
974  center_fiber += (fiber_length / optimal_fiber_length - 1) * 0.5 * vector_fiber;
975  vector_fiber *= fiber_length / optimal_fiber_length;
976 
977  // const G4Vector3D v1_new = center_fiber - 0.5 *vector_fiber;
978 
979  if (_geom->get_construction_verbose() >= 3)
980  cout
981  << "PHG4SpacalPrototype4Detector::Construct_Fibers_SameLengthFiberPerTower::"
982  << GetName() << " - constructed fiber " << fiber_ID << ss.str()
983  //
984  << ", Length = " << optimal_fiber_length << "-"
985  << (optimal_fiber_length - fiber_length) << "mm, " //
986  << "x = " << center_fiber.x() << "mm, " //
987  << "y = " << center_fiber.y() << "mm, " //
988  << "z = " << center_fiber.z() << "mm, " //
989  << "vx = " << vector_fiber.x() << "mm, " //
990  << "vy = " << vector_fiber.y() << "mm, " //
991  << "vz = " << vector_fiber.z() << "mm, " //
992  << endl;
993 
994  const G4double rotation_angle = G4Vector3D(0, 0, 1).angle(vector_fiber);
995  const G4Vector3D rotation_axis =
996  rotation_angle == 0 ? G4Vector3D(1, 0, 0) : G4Vector3D(0, 0, 1).cross(vector_fiber);
997 
998  G4Transform3D fiber_place(
999  G4Translate3D(center_fiber.x(), center_fiber.y(), center_fiber.z()) * G4Rotate3D(rotation_angle, rotation_axis));
1000 
1001  ostringstream name;
1002  name << GetName() + string("_Tower") << g_tower.id << "_fiber"
1003  << ss.str();
1004 
1005  const bool overlapcheck_fiber = OverlapCheck() and (_geom->get_construction_verbose() >= 3);
1006  G4PVPlacement* fiber_physi = new G4PVPlacement(fiber_place, fiber_logic,
1007  G4String(name.str()), LV_tower, false, fiber_ID,
1008  overlapcheck_fiber);
1009  fiber_vol[fiber_physi] = fiber_ID;
1010 
1011  fiber_count++;
1012  }
1013 
1014  if (_geom->get_construction_verbose() >= 2)
1015  cout
1016  << "PHG4SpacalPrototype4Detector::Construct_Fibers_SameLengthFiberPerTower::"
1017  << GetName() << " - constructed tower ID " << g_tower.id << " with "
1018  << fiber_count << " fibers. Average fiber length cut = "
1019  << accumulate(fiber_cut.begin(), fiber_cut.end(), 0.0) / fiber_cut.size() << " mm" << endl;
1020 
1021  return fiber_count;
1022 }
1023 
1025 G4LogicalVolume*
1026 PHG4SpacalPrototype4Detector::Construct_Tower(
1027  const PHG4SpacalPrototype4Detector::SpacalGeom_t::geom_tower& g_tower)
1028 {
1029  assert(_geom);
1030 
1031  if (_geom->get_construction_verbose() >= 2)
1032  {
1033  cout << "PHG4SpacalPrototype4Detector::Construct_Tower::" << GetName()
1034  << " - constructed tower ID " << g_tower.id
1035  << " with geometry parameter: ";
1036  g_tower.identify(cout);
1037  }
1038 
1039  std::ostringstream sout;
1040  sout << "_" << g_tower.id;
1041  const G4String sTowerID(sout.str());
1042 
1043  //Processed PostionSeeds 1 from 1 1
1044 
1045  G4Trap* block_solid = new G4Trap(
1046  /*const G4String& pName*/ G4String(GetName()) + sTowerID,
1047  g_tower.pDz * cm, // G4double pDz,
1048  g_tower.pTheta * rad, g_tower.pPhi * rad, // G4double pTheta, G4double pPhi,
1049  g_tower.pDy1 * cm, g_tower.pDx1 * cm, g_tower.pDx2 * cm, // G4double pDy1, G4double pDx1, G4double pDx2,
1050  g_tower.pAlp1 * rad, // G4double pAlp1,
1051  g_tower.pDy2 * cm, g_tower.pDx3 * cm, g_tower.pDx4 * cm, // G4double pDy2, G4double pDx3, G4double pDx4,
1052  g_tower.pAlp2 * rad // G4double pAlp2 //
1053  );
1054 
1055  G4Material* cylinder_mat = G4Material::GetMaterial(
1056  _geom->get_absorber_mat());
1057  assert(cylinder_mat);
1058 
1059  G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, cylinder_mat,
1060  G4String(G4String(GetName()) + string("_Tower") + sTowerID), 0, 0,
1061  step_limits);
1062 
1063  G4VisAttributes* VisAtt = new G4VisAttributes();
1064  // PHG4Utils::SetColour(VisAtt, "W_Epoxy");
1065  VisAtt->SetColor(.3, .3, .3, .3);
1066  VisAtt->SetVisibility(
1067  _geom->is_azimuthal_seg_visible() or _geom->is_virualize_fiber());
1068  VisAtt->SetForceSolid(not _geom->is_virualize_fiber());
1069  block_logic->SetVisAttributes(VisAtt);
1070 
1071  // construct fibers
1072 
1073  int fiber_count = 0;
1074 
1075  fiber_count = Construct_Fibers_SameLengthFiberPerTower(g_tower, block_logic);
1076 
1077  if (_geom->get_construction_verbose() >= 2)
1078  cout << "PHG4SpacalPrototype4Detector::Construct_Tower::" << GetName()
1079  << " - constructed tower ID " << g_tower.id << " with " << fiber_count
1080  << " fibers using Construct_Fibers_SameLengthFiberPerTower" << endl;
1081 
1082  return block_logic;
1083 }
1084 
1085 G4LogicalVolume*
1086 PHG4SpacalPrototype4Detector::Construct_LightGuide(
1087  const PHG4SpacalPrototype4Detector::SpacalGeom_t::geom_tower& g_tower,
1088  const int index_x, const int index_y)
1089 {
1090  assert(_geom);
1091 
1092  std::ostringstream sout;
1093  sout << "_Lightguide_" << g_tower.id << "_" << index_x << "_" << index_y;
1094  const G4String sTowerID(sout.str());
1095 
1096  assert(g_tower.LightguideHeight > 0);
1097 
1098  // light guide parameters in PHENIX units
1099  const double weight_x1 = 1 - (double) index_y / g_tower.NSubtowerY;
1100  const double weight_x2 = 1 - (double) (index_y + 1) / g_tower.NSubtowerY;
1101  const double weight_xcenter = 1 - (double) (index_y + 0.5) / g_tower.NSubtowerY;
1102 
1103  assert(weight_x1 >= 0 and weight_x1 <= 1);
1104  assert(weight_x2 >= 0 and weight_x2 <= 1);
1105  assert(weight_xcenter >= 0 and weight_xcenter <= 1);
1106 
1107  const double lg_pDx1 = (g_tower.pDx1 * weight_x1 //
1108  + g_tower.pDx2 * (1 - weight_x1)) /
1109  g_tower.NSubtowerX;
1110  const double lg_pDx2 = (g_tower.pDx1 * weight_x2 //
1111  + g_tower.pDx2 * (1 - weight_x2)) /
1112  g_tower.NSubtowerX;
1113  const double lg_pDy1 = g_tower.pDy1 / g_tower.NSubtowerY;
1114  const double lg_Alp1 = atan(
1115  (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));
1116 
1117  const double shift_xcenter = (g_tower.pDx1 * weight_xcenter //
1118  + g_tower.pDx2 * (1 - weight_xcenter)) //
1119  * //
1120  (-g_tower.NSubtowerX + 1. + 2 * index_x) / (double) (g_tower.NSubtowerX);
1121  const double shift_ycenter = g_tower.pDy1 //
1122  * //
1123  (-g_tower.NSubtowerY + 1. + 2 * index_y) / (double) (g_tower.NSubtowerY);
1124 
1125  G4VSolid* block_solid = new G4Trap(
1126  /*const G4String& pName*/ G4String(GetName()) + sTowerID,
1127  0.5 * g_tower.LightguideHeight * cm, // G4double pDz,
1128  0 * rad, 0 * rad, // G4double pTheta, G4double pPhi,
1129  g_tower.LightguideTaperRatio * lg_pDy1 * cm,
1130  g_tower.LightguideTaperRatio * lg_pDx1 * cm,
1131  g_tower.LightguideTaperRatio * lg_pDx2 * cm, // G4double pDy1, G4double pDx1, G4double pDx2,
1132  lg_Alp1 * rad, // G4double pAlp1,
1133  lg_pDy1 * cm, lg_pDx1 * cm, lg_pDx2 * cm, // G4double pDy2, G4double pDx3, G4double pDx4,
1134  lg_Alp1 * rad // G4double pAlp2 //
1135  );
1136 
1137  block_solid = new G4DisplacedSolid(G4String(GetName() + "_displaced"),
1138  block_solid, 0, //
1139  G4ThreeVector( //
1140  tan(g_tower.pTheta * rad) * cos(g_tower.pPhi * rad), //
1141  tan(g_tower.pTheta * rad) * sin(g_tower.pPhi * rad), //
1142  1) * // G4ThreeVector
1143  -(g_tower.pDz) *
1144  cm //
1145  + G4ThreeVector(shift_xcenter * cm, shift_ycenter * cm, 0) // shit in subtower direction
1146  + G4ThreeVector(0, 0, -0.5 * g_tower.LightguideHeight * cm) //shift in the light guide height
1147  );
1148 
1149  G4Material* cylinder_mat = G4Material::GetMaterial(
1150  g_tower.LightguideMaterial);
1151  assert(cylinder_mat);
1152 
1153  G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, cylinder_mat,
1154  G4String(G4String(GetName()) + string("_Tower") + sTowerID), 0, 0,
1155  step_limits);
1156 
1157  G4VisAttributes* VisAtt = new G4VisAttributes();
1158  PHG4Utils::SetColour(VisAtt, g_tower.LightguideMaterial);
1159  // VisAtt->SetColor(.3, .3, .3, .3);
1160  VisAtt->SetVisibility(
1161  _geom->is_azimuthal_seg_visible() or _geom->is_virualize_fiber());
1162  VisAtt->SetForceSolid(not _geom->is_virualize_fiber());
1163  block_logic->SetVisAttributes(VisAtt);
1164 
1165  return block_logic;
1166 }