Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4SpacalPrototypeDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4SpacalPrototypeDetector.cc
1 #include "PHG4SpacalPrototypeDetector.h"
2 
3 #include <g4detectors/PHG4CylinderGeomContainer.h>
4 #include <g4detectors/PHG4CylinderGeom_Spacalv1.h> // for PHG4CylinderGeom_Spacalv1:...
5 
6 #include <phparameter/PHParameters.h>
7 
8 #include <g4detectors/PHG4CylinderGeom_Spacalv3.h> // for PHG4CylinderGeom_...
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/G4DisplacedSolid.hh>
21 #include <Geant4/G4LogicalVolume.hh>
22 #include <Geant4/G4Material.hh>
23 #include <Geant4/G4PVPlacement.hh>
24 #include <Geant4/G4PhysicalConstants.hh>
25 #include <Geant4/G4String.hh> // for G4String
26 #include <Geant4/G4SubtractionSolid.hh>
27 #include <Geant4/G4SystemOfUnits.hh>
28 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
29 #include <Geant4/G4Transform3D.hh> // for G4Transform3D, G4Translate3D
30 #include <Geant4/G4Trap.hh>
31 #include <Geant4/G4Tubs.hh>
32 #include <Geant4/G4Types.hh> // for G4double
33 #include <Geant4/G4UserLimits.hh>
34 #include <Geant4/G4Vector3D.hh>
35 #include <Geant4/G4VisAttributes.hh>
36 
37 #include <boost/foreach.hpp>
38 
39 #include <algorithm>
40 #include <cassert>
41 #include <cmath>
42 #include <iostream> // for operator<<, basic_ostream
43 #include <numeric> // std::accumulate
44 #include <sstream>
45 #include <string> // std::string, std::to_string
46 #include <vector> // for vector
47 
48 class PHG4CylinderGeom;
49 
50 using namespace std;
51 
52 //_______________________________________________________________
53 //note this inactive thickness is ~1.5% of a radiation length
54 PHG4SpacalPrototypeDetector::PHG4SpacalPrototypeDetector(PHG4Subsystem* subsys, PHCompositeNode* Node, PHParameters* parameters, const std::string& dnam)
55  : PHG4Detector(subsys, Node, dnam)
56  , construction_params(parameters)
57  , cylinder_solid(nullptr)
58  , cylinder_logic(nullptr)
59  , cylinder_physi(nullptr)
60  , //
61  active(0)
62  , absorberactive(0)
63  , //
64  step_limits(nullptr)
65  , clading_step_limits(nullptr)
66  , fiber_core_step_limits(nullptr)
67  , //
68  _geom(nullptr)
69 {
70 }
71 
72 PHG4SpacalPrototypeDetector::~PHG4SpacalPrototypeDetector(void)
73 {
74  // deleting nullptr pointers is allowed (results in NOOP)
75  // so checking for not null before deleting is not needed
76  delete step_limits;
77  delete clading_step_limits;
78  delete fiber_core_step_limits;
79  delete _geom;
80 }
81 
82 //_______________________________________________________________
83 int PHG4SpacalPrototypeDetector::IsInCylinderActive(
84  const G4VPhysicalVolume* volume)
85 {
86  // cout << "checking detector" << endl;
87  if (active && fiber_core_vol.find(volume) != fiber_core_vol.end())
88  {
89  // return fiber_core_vol.find(volume)->second;
90  return FIBER_CORE;
91  }
92  if (absorberactive)
93  {
94  if (fiber_vol.find(volume) != fiber_vol.end())
95  return FIBER_CLADING;
96 
97  if (block_vol.find(volume) != block_vol.end())
98  return ABSORBER;
99 
100  if (calo_vol.find(volume) != calo_vol.end())
101  return SUPPORT;
102  }
103  return INACTIVE;
104 }
105 
106 //_______________________________________________________________
107 void PHG4SpacalPrototypeDetector::ConstructMe(G4LogicalVolume* logicWorld)
108 {
109  PHNodeIterator iter(topNode());
110  PHCompositeNode* parNode = dynamic_cast<PHCompositeNode*>(iter.findFirst(
111  "PHCompositeNode", "RUN"));
112  assert(parNode);
113 
114  if (!_geom)
115  {
116  _geom = new SpacalGeom_t();
117  }
118  assert(_geom);
119 
120  _geom->set_config(
121  SpacalGeom_t::kFullProjective_2DTaper_SameLengthFiberPerTower);
122  // _geom->load_demo_sector_tower_map4();
123  // _geom->set_construction_verbose(2);
124  _geom->ImportParameters(*construction_params);
125 
126  _geom->subtower_consistency_check();
127 
128  // step_limits = new G4UserLimits(_geom->get_calo_step_size() * cm);
129  //
130  // clading_step_limits = new G4UserLimits(
131  // _geom->get_fiber_clading_step_size() * cm);
132  step_limits = nullptr;
133  clading_step_limits = nullptr;
134 
135  fiber_core_step_limits = new G4UserLimits(
136  _geom->get_fiber_core_step_size() * cm);
137 
138  const G4double z_rotation = construction_params->get_double_param(
139  "z_rotation_degree") *
140  degree;
141  const G4double enclosure_x = construction_params->get_double_param(
142  "enclosure_x") *
143  cm;
144  const G4double enclosure_x_shift = construction_params->get_double_param(
145  "enclosure_x_shift") *
146  cm;
147  const G4double enclosure_y = construction_params->get_double_param(
148  "enclosure_y") *
149  cm;
150  const G4double enclosure_z = construction_params->get_double_param(
151  "enclosure_z") *
152  cm;
153 
154  const G4double box_x_shift = (_geom->get_radius() + 0.5 * _geom->get_thickness()) * cm + enclosure_x_shift;
155  const G4double box_z_shift = 0.5 * (_geom->get_zmin() + _geom->get_zmax()) * cm;
156 
157  // G4Tubs* _cylinder_solid = new G4Tubs(G4String(GetName()),
158  // _geom->get_radius() * cm, _geom->get_max_radius() * cm,
159  // _geom->get_length() * cm / 2.0, 0, twopi);
160  // G4VSolid* cylinder_solid = new G4Box(G4String(GetName()),
161  // _geom->get_thickness() * 1.1 * 0.5 * cm, //
162  // twopi / _geom->get_azimuthal_n_sec() * _geom->get_sector_map().size()
163  // * 0.5 * _geom->get_max_radius() * cm, //
164  // _geom->get_length() * cm / 2.0);
165  G4VSolid* cylinder_solid = new G4Box(G4String(GetName()),
166  enclosure_x * 0.5, //
167  enclosure_y * 0.5, //
168  enclosure_z * 0.5);
169 
170  cylinder_solid = new G4DisplacedSolid(G4String(GetName() + "_displaced"),
171  cylinder_solid, 0, //
172  G4ThreeVector(box_x_shift, 0, box_z_shift));
173 
174  G4Material* cylinder_mat = G4Material::GetMaterial("G4_AIR");
175  assert(cylinder_mat);
176 
177  cylinder_logic = new G4LogicalVolume(cylinder_solid, cylinder_mat,
178  G4String(GetName()), 0, 0, 0);
179  G4VisAttributes* VisAtt = new G4VisAttributes();
180  PHG4Utils::SetColour(VisAtt, "W_Epoxy");
181  VisAtt->SetVisibility(true);
182  VisAtt->SetForceSolid(
183  (not _geom->is_virualize_fiber()) and (not _geom->is_azimuthal_seg_visible()));
184  cylinder_logic->SetVisAttributes(VisAtt);
185 
186  G4Transform3D cylinder_place(
187  G4Translate3D(_geom->get_xpos() * cm, _geom->get_ypos() * cm,
188  _geom->get_zpos() * cm) //
189  * G4RotateZ3D(z_rotation) //
190  * G4Translate3D(
191  -(_geom->get_radius() + 0.5 * _geom->get_thickness()) * cm, 0,
192  0));
193 
194  cylinder_physi = new G4PVPlacement(cylinder_place, cylinder_logic,
195  G4String(GetName()), logicWorld, false, 0, OverlapCheck());
196 
197  // install sectors
198  std::pair<G4LogicalVolume*, G4Transform3D> psec = Construct_AzimuthalSeg();
199  G4LogicalVolume* sec_logic = psec.first;
200  const G4Transform3D& sec_trans = psec.second;
201  BOOST_FOREACH (const SpacalGeom_t::sector_map_t::value_type& val, _geom->get_sector_map())
202  {
203  const int sec = val.first;
204  const double rot = val.second;
205 
206  G4Transform3D sec_place = G4RotateZ3D(rot) * sec_trans;
207 
208  ostringstream name;
209  name << GetName() << "_sec" << sec;
210 
211  G4PVPlacement* calo_phys = new G4PVPlacement(sec_place, sec_logic,
212  G4String(name.str()), cylinder_logic, false, sec,
213  OverlapCheck());
214  calo_vol[calo_phys] = sec;
215  }
216  _geom->set_nscint(_geom->get_nscint() * _geom->get_sector_map().size());
217 
218  // install electronics
219  const G4double electronics_thickness = construction_params->get_double_param(
220  "electronics_thickness") *
221  cm;
222  const string electronics_material = construction_params->get_string_param(
223  "electronics_material");
224 
225  G4VSolid* electronics_solid = new G4Box(G4String(GetName()),
226  electronics_thickness / 2.0, //
227  sin(
228  twopi / _geom->get_azimuthal_n_sec() * _geom->get_sector_map().size() / 2) *
229  (_geom->get_radius() - _geom->get_max_lightguide_height()) * cm, //
230  _geom->get_length() * cm / 2.0);
231 
232  electronics_solid = new G4DisplacedSolid(G4String(GetName() + "_displaced"),
233  electronics_solid,
234  0, //
235  G4ThreeVector(
236  cos(
237  twopi / _geom->get_azimuthal_n_sec() * _geom->get_sector_map().size() / 2) *
238  (_geom->get_radius() - _geom->get_max_lightguide_height()) * cm -
239  electronics_thickness,
240  0, box_z_shift));
241 
242  G4Material* electronics_mat = G4Material::GetMaterial(electronics_material);
243  assert(electronics_mat);
244 
245  G4LogicalVolume* electronics_logic = new G4LogicalVolume(electronics_solid,
246  electronics_mat, G4String(GetName()) + "_Electronics", 0, 0, 0);
247  VisAtt = new G4VisAttributes();
248  PHG4Utils::SetColour(VisAtt, electronics_material);
249  VisAtt->SetVisibility(true);
250  VisAtt->SetForceSolid(not _geom->is_virualize_fiber());
251  electronics_logic->SetVisAttributes(VisAtt);
252 
253  new G4PVPlacement(G4Transform3D::Identity, electronics_logic,
254  G4String(GetName()) + "_Electronics", cylinder_logic, false, 0,
255  OverlapCheck());
256 
257  // install the enclosure box
258  const G4double enclosure_thickness = construction_params->get_double_param(
259  "enclosure_thickness") *
260  cm;
261  const string enclosure_material = construction_params->get_string_param(
262  "enclosure_material");
263 
264  G4VSolid* outer_enclosur_solid = new G4Box(
265  G4String(GetName()) + "_outer_enclosur_solid", enclosure_x * 0.5, //
266  enclosure_y * 0.5, //
267  enclosure_z * 0.5);
268  G4VSolid* inner_enclosur_solid = new G4Box(
269  G4String(GetName()) + "_inner_enclosur_solid",
270  enclosure_x * 0.5 - enclosure_thickness, //
271  enclosure_y * 0.5 - enclosure_thickness, //
272  enclosure_z * 0.5 - enclosure_thickness);
273  G4VSolid* enclosure_solid = new G4SubtractionSolid(
274  G4String(GetName()) + "_enclosure_solid", //
275  outer_enclosur_solid, inner_enclosur_solid);
276  enclosure_solid = new G4DisplacedSolid(
277  G4String(GetName() + "_enclosure_solid_displaced"), enclosure_solid, 0, //
278  G4ThreeVector(box_x_shift, 0, box_z_shift));
279 
280  G4Material* enclosure_mat = G4Material::GetMaterial(enclosure_material);
281  assert(enclosure_mat);
282 
283  G4LogicalVolume* enclosure_logic = new G4LogicalVolume(enclosure_solid,
284  enclosure_mat, G4String(GetName()) + "_enclosure", 0, 0, 0);
285  VisAtt = new G4VisAttributes();
286  PHG4Utils::SetColour(VisAtt, enclosure_material);
287  VisAtt->SetVisibility(true);
288  VisAtt->SetForceSolid(not _geom->is_virualize_fiber());
289  enclosure_logic->SetVisAttributes(VisAtt);
290 
291  new G4PVPlacement(G4Transform3D::Identity, enclosure_logic,
292  G4String(GetName()) + "_enclosure", cylinder_logic, false, 0,
293  OverlapCheck());
294 
295  if (active)
296  {
297  ostringstream geonode;
298  if (superdetector != "NONE")
299  {
300  geonode << "CYLINDERGEOM_" << superdetector;
301  }
302  else
303  {
304  geonode << "CYLINDERGEOM_" << detector_type;
305  }
306  PHG4CylinderGeomContainer* geo = findNode::getClass<
307  PHG4CylinderGeomContainer>(topNode(), geonode.str());
308  if (!geo)
309  {
310  geo = new PHG4CylinderGeomContainer();
311  PHNodeIterator iter(topNode());
312  PHCompositeNode* runNode =
313  dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode",
314  "RUN"));
315  PHIODataNode<PHObject>* newNode = new PHIODataNode<PHObject>(geo,
316  geonode.str(), "PHObject");
317  runNode->addNode(newNode);
318  }
319  // here in the detector class we have internal units, convert to cm
320  // before putting into the geom object
321  PHG4CylinderGeom* mygeom = new SpacalGeom_t(*_geom);
322  geo->AddLayerGeom(0, mygeom);
323  if (_geom->get_construction_verbose() >= 1)
324  {
325  cout << "PHG4SpacalPrototypeDetector::Construct::" << GetName()
326  << " - Print Layer Geometry:" << endl;
327  geo->identify();
328  }
329  }
330 
331  if (absorberactive)
332  {
333  ostringstream geonode;
334  if (superdetector != "NONE")
335  {
336  geonode << "CYLINDERGEOM_ABSORBER_" << superdetector;
337  }
338  else
339  {
340  geonode << "CYLINDERGEOM_ABSORBER_" << detector_type << "_" << 0;
341  }
342  PHG4CylinderGeomContainer* geo = findNode::getClass<PHG4CylinderGeomContainer>(topNode(), geonode.str());
343  if (!geo)
344  {
345  geo = new PHG4CylinderGeomContainer();
346  PHNodeIterator iter(topNode());
347  PHCompositeNode* runNode =
348  dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode",
349  "RUN"));
350  PHIODataNode<PHObject>* newNode = new PHIODataNode<PHObject>(geo, geonode.str(), "PHObject");
351  runNode->addNode(newNode);
352  }
353  // here in the detector class we have internal units, convert to cm
354  // before putting into the geom object
355  PHG4CylinderGeom* mygeom = new SpacalGeom_t(*_geom);
356  geo->AddLayerGeom(0, mygeom);
357  if (_geom->get_construction_verbose() >= 1)
358  {
359  cout << "PHG4SpacalPrototypeDetector::Construct::" << GetName()
360  << " - Print Absorber Layer Geometry:" << endl;
361  geo->identify();
362  }
363  }
364 
365  if (_geom->get_construction_verbose() >= 1)
366  {
367  cout << "PHG4SpacalPrototypeDetector::Construct::" << GetName()
368  << " - Completed. Print G4 Geometry:" << endl;
369  Print();
370  cout << "box_x_shift = " << box_x_shift << endl;
371  cout << "box_z_shift = " << box_z_shift << endl;
372  }
373 }
374 
375 std::pair<G4LogicalVolume*, G4Transform3D>
376 PHG4SpacalPrototypeDetector::Construct_AzimuthalSeg()
377 {
378  assert(_geom);
379  assert(_geom->get_azimuthal_n_sec() > 4);
380 
381  G4VSolid* sec_solid = new G4Tubs(G4String(GetName() + string("_sec")),
382  (_geom->get_radius() - _geom->get_max_lightguide_height()) * cm,
383  _geom->get_max_radius() * cm, _geom->get_length() * cm / 2.0,
384  halfpi - pi / _geom->get_azimuthal_n_sec(),
385  twopi / _geom->get_azimuthal_n_sec());
386 
387  const G4double box_z_shift = 0.5 * (_geom->get_zmin() + _geom->get_zmax()) * cm;
388  sec_solid = new G4DisplacedSolid(G4String(GetName() + "_displaced" + string("_sec")),
389  sec_solid, 0, //
390  G4ThreeVector(0, 0, box_z_shift));
391 
392  G4Material* cylinder_mat = G4Material::GetMaterial("G4_AIR");
393  assert(cylinder_mat);
394 
395  G4LogicalVolume* sec_logic = new G4LogicalVolume(sec_solid, cylinder_mat,
396  G4String(G4String(GetName() + string("_sec"))), 0, 0, step_limits);
397 
398  G4VisAttributes* VisAtt = new G4VisAttributes();
399  VisAtt->SetColor(.5, .9, .5, .1);
400  VisAtt->SetVisibility(
401  _geom->is_azimuthal_seg_visible() or _geom->is_virualize_fiber());
402  VisAtt->SetForceSolid(false);
403  VisAtt->SetForceWireframe(true);
404  sec_logic->SetVisAttributes(VisAtt);
405 
406  // // construct walls
407  // G4Material * wall_mat = G4Material::GetMaterial(_geom->get_sidewall_mat());
408  // assert(wall_mat);
409  //
410  // G4VisAttributes* wall_VisAtt = new G4VisAttributes();
411  // VisAtt->SetColor(.5, .9, .5, .5);
412  // wall_VisAtt->SetVisibility(
413  // _geom->is_azimuthal_seg_visible() and (not _geom->is_virualize_fiber()));
414  // wall_VisAtt->SetForceSolid(true);
415  //
416  assert(_geom->get_sidewall_thickness() == 0);
417  //section skin not used in prototype construction
418  // if (_geom->get_sidewall_thickness() > 0)
419  // {
420  // // end walls
421  // if (_geom->get_construction_verbose() >= 1)
422  // {
423  // cout << "PHG4SpacalPrototypeDetector::Construct_AzimuthalSeg::"
424  // << GetName() << " - construct end walls." << endl;
425  // }
426  // G4Tubs* wall_solid = new G4Tubs(G4String(GetName() + string("_EndWall")),
427  // _geom->get_radius() * cm + _geom->get_sidewall_outer_torr() * cm,
428  // _geom->get_max_radius() * cm - _geom->get_sidewall_outer_torr() * cm,
429  // _geom->get_sidewall_thickness() * cm / 2.0,
430  // halfpi - pi / _geom->get_azimuthal_n_sec(),
431  // twopi / _geom->get_azimuthal_n_sec());
432  //
433  // G4LogicalVolume * wall_logic = new G4LogicalVolume(wall_solid, wall_mat,
434  // G4String(G4String(GetName() + string("_EndWall"))), 0, 0,
435  // step_limits);
436  // wall_logic->SetVisAttributes(wall_VisAtt);
437  //
438  // typedef map<int, double> z_locations_t;
439  // z_locations_t z_locations;
440  // z_locations[1000] = _geom->get_sidewall_thickness() * cm / 2.0
441  // + _geom->get_assembly_spacing() * cm;
442  // z_locations[1001] = _geom->get_length() * cm / 2.0
443  // - (_geom->get_sidewall_thickness() * cm / 2.0
444  // + _geom->get_assembly_spacing() * cm);
445  // z_locations[1100] = -(_geom->get_sidewall_thickness() * cm / 2.0
446  // + _geom->get_assembly_spacing() * cm);
447  // z_locations[1101] = -(_geom->get_length() * cm / 2.0
448  // - (_geom->get_sidewall_thickness() * cm / 2.0
449  // + _geom->get_assembly_spacing() * cm));
450  //
451  // BOOST_FOREACH(z_locations_t::value_type& val, z_locations)
452  // {
453  // if (_geom->get_construction_verbose() >= 2)
454  // cout << "PHG4SpacalPrototypeDetector::Construct_AzimuthalSeg::"
455  // << GetName() << " - constructed End Wall ID " << val.first
456  // << " @ Z = " << val.second << endl;
457  //
458  // G4Transform3D wall_trans = G4TranslateZ3D(val.second);
459  //
460  // G4PVPlacement * wall_phys = new G4PVPlacement(wall_trans, wall_logic,
461  // G4String(GetName()) + G4String("_EndWall"), sec_logic,
462  // false, val.first, OverlapCheck());
463  //
464  // calo_vol[wall_phys] = val.first;
465  // }
466  // }
467  //
468  // if (_geom->get_sidewall_thickness() > 0)
469  // {
470  // // side walls
471  // if (_geom->get_construction_verbose() >= 1)
472  // {
473  // cout << "PHG4SpacalPrototypeDetector::Construct_AzimuthalSeg::"
474  // << GetName() << " - construct side walls." << endl;
475  // }
476  // G4Box* wall_solid = new G4Box(G4String(GetName() + string("_SideWall")),
477  // _geom->get_sidewall_thickness() * cm / 2.0,
478  // _geom->get_thickness() * cm / 2.
479  // - 2 * _geom->get_sidewall_outer_torr() * cm,
480  // (_geom->get_length() / 2.
481  // - 2
482  // * (_geom->get_sidewall_thickness()
483  // + 2. * _geom->get_assembly_spacing())) * cm * .5);
484  //
485  // G4LogicalVolume * wall_logic = new G4LogicalVolume(wall_solid, wall_mat,
486  // G4String(G4String(GetName() + string("_SideWall"))), 0, 0,
487  // step_limits);
488  // wall_logic->SetVisAttributes(wall_VisAtt);
489  //
490  // typedef map<int, pair<int, int> > sign_t;
491  // sign_t signs;
492  // signs[2000] = make_pair(+1, +1);
493  // signs[2001] = make_pair(+1, -1);
494  // signs[2100] = make_pair(-1, +1);
495  // signs[2101] = make_pair(-1, -1);
496  //
497  // BOOST_FOREACH(sign_t::value_type& val, signs)
498  // {
499  // const int sign_z = val.second.first;
500  // const int sign_azimuth = val.second.second;
501  //
502  // if (_geom->get_construction_verbose() >= 2)
503  // cout << "PHG4SpacalPrototypeDetector::Construct_AzimuthalSeg::"
504  // << GetName() << " - constructed Side Wall ID " << val.first
505  // << " with" << " Shift X = "
506  // << sign_azimuth
507  // * (_geom->get_sidewall_thickness() * cm / 2.0
508  // + _geom->get_sidewall_outer_torr() * cm)
509  // << " Rotation Z = "
510  // << sign_azimuth * pi / _geom->get_azimuthal_n_sec()
511  // << " Shift Z = " << sign_z * (_geom->get_length() * cm / 4)
512  // << endl;
513  //
514  // G4Transform3D wall_trans = G4RotateZ3D(
515  // sign_azimuth * pi / _geom->get_azimuthal_n_sec())
516  // * G4TranslateZ3D(sign_z * (_geom->get_length() * cm / 4))
517  // * G4TranslateY3D(
518  // _geom->get_radius() * cm + _geom->get_thickness() * cm / 2.)
519  // * G4TranslateX3D(
520  // sign_azimuth
521  // * (_geom->get_sidewall_thickness() * cm / 2.0
522  // + _geom->get_sidewall_outer_torr() * cm));
523  //
524  // G4PVPlacement * wall_phys = new G4PVPlacement(wall_trans, wall_logic,
525  // G4String(GetName()) + G4String("_EndWall"), sec_logic,
526  // false, val.first, OverlapCheck());
527  //
528  // calo_vol[wall_phys] = val.first;
529  // }
530  // }
531 
532  // construct towers
533 
534  BOOST_FOREACH (const SpacalGeom_t::tower_map_t::value_type& val, _geom->get_sector_tower_map())
535  {
536  const SpacalGeom_t::geom_tower& g_tower = val.second;
537  G4LogicalVolume* LV_tower = Construct_Tower(g_tower);
538 
539  G4Transform3D block_trans = G4TranslateX3D(g_tower.centralX * cm) * G4TranslateY3D(g_tower.centralY * cm) * G4TranslateZ3D(g_tower.centralZ * cm) * G4RotateX3D(g_tower.pRotationAngleX * rad);
540 
541  const bool overlapcheck_block = OverlapCheck() and (_geom->get_construction_verbose() >= 2);
542 
543  G4PVPlacement* block_phys = new G4PVPlacement(block_trans, LV_tower,
544  LV_tower->GetName(), sec_logic, false, g_tower.id,
545  overlapcheck_block);
546  block_vol[block_phys] = g_tower.id;
547 
548  if (g_tower.LightguideHeight > 0)
549  {
550  // also build a light guide
551 
552  for (int ix = 0; ix < g_tower.NSubtowerX; ix++)
553  // int ix = 0;
554  {
555  for (int iy = 0; iy < g_tower.NSubtowerY; iy++)
556  // int iy = 0;
557  {
558  G4LogicalVolume* LV_lg = Construct_LightGuide(g_tower, ix,
559  iy);
560 
561  new G4PVPlacement(block_trans, LV_lg, LV_lg->GetName(),
562  sec_logic, false, g_tower.id, overlapcheck_block);
563  }
564  }
565  }
566  }
567 
568  cout << "PHG4SpacalPrototypeDetector::Construct_AzimuthalSeg::" << GetName()
569  << " - constructed " << _geom->get_sector_tower_map().size()
570  << " unique towers" << endl;
571 
572  return make_pair(sec_logic, G4Transform3D::Identity);
573 }
574 
575 G4LogicalVolume*
576 PHG4SpacalPrototypeDetector::Construct_Fiber(const G4double length,
577  const string& id)
578 {
579  G4Tubs* fiber_solid = new G4Tubs(G4String(GetName() + string("_fiber") + id),
580  0, _geom->get_fiber_outer_r() * cm, length / 2.0, 0, twopi);
581 
582  G4Material* clading_mat = G4Material::GetMaterial(
583  _geom->get_fiber_clading_mat());
584  assert(clading_mat);
585 
586  G4LogicalVolume* fiber_logic = new G4LogicalVolume(fiber_solid, clading_mat,
587  G4String(G4String(GetName() + string("_fiber") + id)), 0, 0,
588  clading_step_limits);
589 
590  {
591  G4VisAttributes* VisAtt = new G4VisAttributes();
592  PHG4Utils::SetColour(VisAtt, "G4_POLYSTYRENE");
593  VisAtt->SetVisibility(_geom->is_virualize_fiber());
594  VisAtt->SetForceSolid(_geom->is_virualize_fiber());
595  fiber_logic->SetVisAttributes(VisAtt);
596  }
597 
598  G4Tubs* core_solid = new G4Tubs(
599  G4String(GetName() + string("_fiber_core") + id), 0,
600  _geom->get_fiber_core_diameter() * cm / 2, length / 2.0, 0, twopi);
601 
602  G4Material* core_mat = G4Material::GetMaterial(_geom->get_fiber_core_mat());
603  assert(core_mat);
604 
605  G4LogicalVolume* core_logic = new G4LogicalVolume(core_solid, core_mat,
606  G4String(G4String(GetName() + string("_fiber_core") + id)), 0, 0,
607  fiber_core_step_limits);
608 
609  {
610  G4VisAttributes* VisAtt = new G4VisAttributes();
611  PHG4Utils::SetColour(VisAtt, "G4_POLYSTYRENE");
612  VisAtt->SetVisibility(false);
613  VisAtt->SetForceSolid(false);
614  core_logic->SetVisAttributes(VisAtt);
615  }
616 
617  const bool overlapcheck_fiber = OverlapCheck() and (_geom->get_construction_verbose() >= 3);
618  G4PVPlacement* core_physi = new G4PVPlacement(0, G4ThreeVector(), core_logic,
619  G4String(G4String(GetName() + string("_fiber_core") + id)), fiber_logic,
620  false, 0, overlapcheck_fiber);
621  fiber_core_vol[core_physi] = 0;
622 
623  return fiber_logic;
624 }
625 
626 void PHG4SpacalPrototypeDetector::Print(const std::string& /*what*/) const
627 {
628  assert(_geom);
629 
630  cout << "PHG4SpacalPrototypeDetector::Print::" << GetName()
631  << " - Print Geometry:" << endl;
632  _geom->Print();
633 
634  return;
635 }
636 
638 int PHG4SpacalPrototypeDetector::Construct_Fibers_SameLengthFiberPerTower(
639  const PHG4SpacalPrototypeDetector::SpacalGeom_t::geom_tower& g_tower,
640  G4LogicalVolume* LV_tower)
641 {
642  assert(_geom);
643 
644  // construct fibers
645 
646  // first check out the fibers geometry
647 
648  typedef map<int, pair<G4Vector3D, G4Vector3D> > fiber_par_map;
649  fiber_par_map fiber_par;
650  G4double min_fiber_length = g_tower.pDz * cm * 4;
651 
652  G4Vector3D v_zshift = G4Vector3D(tan(g_tower.pTheta) * cos(g_tower.pPhi),
653  tan(g_tower.pTheta) * sin(g_tower.pPhi), 1) *
654  g_tower.pDz;
655  // int fiber_ID = 0;
656  for (int ix = 0; ix < g_tower.NFiberX; ix++)
657  // int ix = 0;
658  {
659  const double weighted_ix = static_cast<double>(ix) / (g_tower.NFiberX - 1.);
660 
661  const double weighted_pDx1 = (g_tower.pDx1 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
662  const double weighted_pDx2 = (g_tower.pDx2 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
663 
664  const double weighted_pDx3 = (g_tower.pDx3 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
665  const double weighted_pDx4 = (g_tower.pDx4 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
666 
667  for (int iy = 0; iy < g_tower.NFiberY; iy++)
668  // int iy = 0;
669  {
670  if ((ix + iy) % 2 == 1)
671  continue; // make a triangle pattern
672 
673  const double weighted_iy = static_cast<double>(iy) / (g_tower.NFiberY - 1.);
674 
675  const double weighted_pDy1 = (g_tower.pDy1 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
676  const double weighted_pDy2 = (g_tower.pDy2 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
677 
678  const double weighted_pDx12 = weighted_pDx1 * (1 - weighted_iy) + weighted_pDx2 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp1);
679  const double weighted_pDx34 = weighted_pDx3 * (1 - weighted_iy) + weighted_pDx4 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp2);
680 
681  G4Vector3D v1 = G4Vector3D(weighted_pDx12, weighted_pDy1, 0) - v_zshift;
682  G4Vector3D v2 = G4Vector3D(weighted_pDx34, weighted_pDy2, 0) + v_zshift;
683 
684  G4Vector3D vector_fiber = (v2 - v1);
685  vector_fiber *= (vector_fiber.mag() - _geom->get_fiber_outer_r()) / vector_fiber.mag(); // shrink by fiber boundary protection
686  G4Vector3D center_fiber = (v2 + v1) / 2;
687 
688  // convert to Geant4 units
689  vector_fiber *= cm;
690  center_fiber *= cm;
691 
692  const int fiber_ID = g_tower.compose_fiber_id(ix, iy);
693  fiber_par[fiber_ID] = make_pair(vector_fiber, center_fiber);
694 
695  const G4double fiber_length = vector_fiber.mag();
696 
697  min_fiber_length = min(fiber_length, min_fiber_length);
698 
699  // ++fiber_ID;
700  }
701  }
702 
703  int fiber_count = 0;
704 
705  const G4double fiber_length = min_fiber_length;
706  vector<G4double> fiber_cut;
707 
708  ostringstream ss;
709  ss << string("_Tower") << g_tower.id;
710  G4LogicalVolume* fiber_logic = Construct_Fiber(fiber_length, ss.str());
711 
712  BOOST_FOREACH (const fiber_par_map::value_type& val, fiber_par)
713  {
714  const int fiber_ID = val.first;
715  G4Vector3D vector_fiber = val.second.first;
716  G4Vector3D center_fiber = val.second.second;
717  const G4double optimal_fiber_length = vector_fiber.mag();
718 
719  const G4Vector3D v1 = center_fiber - 0.5 * vector_fiber;
720 
721  // keep a statistics
722  assert(optimal_fiber_length - fiber_length >= 0);
723  fiber_cut.push_back(optimal_fiber_length - fiber_length);
724 
725  center_fiber += (fiber_length / optimal_fiber_length - 1) * 0.5 * vector_fiber;
726  vector_fiber *= fiber_length / optimal_fiber_length;
727 
728  // const G4Vector3D v1_new = center_fiber - 0.5 *vector_fiber;
729 
730  if (_geom->get_construction_verbose() >= 3)
731  cout
732  << "PHG4SpacalPrototypeDetector::Construct_Fibers_SameLengthFiberPerTower::"
733  << GetName() << " - constructed fiber " << fiber_ID << ss.str()
734  //
735  << ", Length = " << optimal_fiber_length << "-"
736  << (optimal_fiber_length - fiber_length) << "mm, " //
737  << "x = " << center_fiber.x() << "mm, " //
738  << "y = " << center_fiber.y() << "mm, " //
739  << "z = " << center_fiber.z() << "mm, " //
740  << "vx = " << vector_fiber.x() << "mm, " //
741  << "vy = " << vector_fiber.y() << "mm, " //
742  << "vz = " << vector_fiber.z() << "mm, " //
743  << endl;
744 
745  const G4double rotation_angle = G4Vector3D(0, 0, 1).angle(vector_fiber);
746  const G4Vector3D rotation_axis =
747  rotation_angle == 0 ? G4Vector3D(1, 0, 0) : G4Vector3D(0, 0, 1).cross(vector_fiber);
748 
749  G4Transform3D fiber_place(
750  G4Translate3D(center_fiber.x(), center_fiber.y(), center_fiber.z()) * G4Rotate3D(rotation_angle, rotation_axis));
751 
752  ostringstream name;
753  name << GetName() + string("_Tower") << g_tower.id << "_fiber"
754  << ss.str();
755 
756  const bool overlapcheck_fiber = OverlapCheck() and (_geom->get_construction_verbose() >= 3);
757  G4PVPlacement* fiber_physi = new G4PVPlacement(fiber_place, fiber_logic,
758  G4String(name.str()), LV_tower, false, fiber_ID,
759  overlapcheck_fiber);
760  fiber_vol[fiber_physi] = fiber_ID;
761 
762  fiber_count++;
763  }
764 
765  if (_geom->get_construction_verbose() >= 2)
766  cout
767  << "PHG4SpacalPrototypeDetector::Construct_Fibers_SameLengthFiberPerTower::"
768  << GetName() << " - constructed tower ID " << g_tower.id << " with "
769  << fiber_count << " fibers. Average fiber length cut = "
770  << accumulate(fiber_cut.begin(), fiber_cut.end(), 0.0) / fiber_cut.size() << " mm" << endl;
771 
772  return fiber_count;
773 }
774 
776 G4LogicalVolume*
777 PHG4SpacalPrototypeDetector::Construct_Tower(
778  const PHG4SpacalPrototypeDetector::SpacalGeom_t::geom_tower& g_tower)
779 {
780  assert(_geom);
781 
782  if (_geom->get_construction_verbose() >= 2)
783  {
784  cout << "PHG4SpacalPrototypeDetector::Construct_Tower::" << GetName()
785  << " - constructed tower ID " << g_tower.id
786  << " with geometry parameter: ";
787  g_tower.identify(cout);
788  }
789 
790  std::ostringstream sout;
791  sout << "_" << g_tower.id;
792  const G4String sTowerID(sout.str());
793 
794  //Processed PostionSeeds 1 from 1 1
795 
796  G4Trap* block_solid = new G4Trap(
797  /*const G4String& pName*/ G4String(GetName()) + sTowerID,
798  g_tower.pDz * cm, // G4double pDz,
799  g_tower.pTheta * rad, g_tower.pPhi * rad, // G4double pTheta, G4double pPhi,
800  g_tower.pDy1 * cm, g_tower.pDx1 * cm, g_tower.pDx2 * cm, // G4double pDy1, G4double pDx1, G4double pDx2,
801  g_tower.pAlp1 * rad, // G4double pAlp1,
802  g_tower.pDy2 * cm, g_tower.pDx3 * cm, g_tower.pDx4 * cm, // G4double pDy2, G4double pDx3, G4double pDx4,
803  g_tower.pAlp2 * rad // G4double pAlp2 //
804  );
805 
806  G4Material* cylinder_mat = G4Material::GetMaterial(
807  _geom->get_absorber_mat());
808  assert(cylinder_mat);
809 
810  G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, cylinder_mat,
811  G4String(G4String(GetName()) + string("_Tower") + sTowerID), 0, 0,
812  step_limits);
813 
814  G4VisAttributes* VisAtt = new G4VisAttributes();
815  // PHG4Utils::SetColour(VisAtt, "W_Epoxy");
816  VisAtt->SetColor(.3, .3, .3, .3);
817  VisAtt->SetVisibility(
818  _geom->is_azimuthal_seg_visible() or _geom->is_virualize_fiber());
819  VisAtt->SetForceSolid(not _geom->is_virualize_fiber());
820  block_logic->SetVisAttributes(VisAtt);
821 
822  // construct fibers
823 
824  int fiber_count = 0;
825 
826  fiber_count = Construct_Fibers_SameLengthFiberPerTower(g_tower, block_logic);
827 
828  if (_geom->get_construction_verbose() >= 2)
829  cout << "PHG4SpacalPrototypeDetector::Construct_Tower::" << GetName()
830  << " - constructed tower ID " << g_tower.id << " with " << fiber_count
831  << " fibers using Construct_Fibers_SameLengthFiberPerTower" << endl;
832 
833  return block_logic;
834 }
835 
836 G4LogicalVolume*
837 PHG4SpacalPrototypeDetector::Construct_LightGuide(
838  const PHG4SpacalPrototypeDetector::SpacalGeom_t::geom_tower& g_tower,
839  const int index_x, const int index_y)
840 {
841  assert(_geom);
842 
843  std::ostringstream sout;
844  sout << "_Lightguide_" << g_tower.id << "_" << index_x << "_" << index_y;
845  const G4String sTowerID(sout.str());
846 
847  assert(g_tower.LightguideHeight > 0);
848 
849  // light guide parameters in PHENIX units
850  const double weight_x1 = 1 - (double) index_y / g_tower.NSubtowerY;
851  const double weight_x2 = 1 - (double) (index_y + 1) / g_tower.NSubtowerY;
852  const double weight_xcenter = 1 - (double) (index_y + 0.5) / g_tower.NSubtowerY;
853 
854  assert(weight_x1 >= 0 and weight_x1 <= 1);
855  assert(weight_x2 >= 0 and weight_x2 <= 1);
856  assert(weight_xcenter >= 0 and weight_xcenter <= 1);
857 
858  const double lg_pDx1 = (g_tower.pDx1 * weight_x1 //
859  + g_tower.pDx2 * (1 - weight_x1)) /
860  g_tower.NSubtowerX;
861  const double lg_pDx2 = (g_tower.pDx1 * weight_x2 //
862  + g_tower.pDx2 * (1 - weight_x2)) /
863  g_tower.NSubtowerX;
864  const double lg_pDy1 = g_tower.pDy1 / g_tower.NSubtowerY;
865  const double lg_Alp1 = atan(
866  (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));
867 
868  const double shift_xcenter = (g_tower.pDx1 * weight_xcenter //
869  + g_tower.pDx2 * (1 - weight_xcenter)) //
870  * //
871  (-g_tower.NSubtowerX + 1. + 2 * index_x) / (double) (g_tower.NSubtowerX);
872  const double shift_ycenter = g_tower.pDy1 //
873  * //
874  (-g_tower.NSubtowerY + 1. + 2 * index_y) / (double) (g_tower.NSubtowerY);
875 
876  G4VSolid* block_solid = new G4Trap(
877  /*const G4String& pName*/ G4String(GetName()) + sTowerID,
878  0.5 * g_tower.LightguideHeight * cm, // G4double pDz,
879  0 * rad, 0 * rad, // G4double pTheta, G4double pPhi,
880  g_tower.LightguideTaperRatio * lg_pDy1 * cm,
881  g_tower.LightguideTaperRatio * lg_pDx1 * cm,
882  g_tower.LightguideTaperRatio * lg_pDx2 * cm, // G4double pDy1, G4double pDx1, G4double pDx2,
883  lg_Alp1 * rad, // G4double pAlp1,
884  lg_pDy1 * cm, lg_pDx1 * cm, lg_pDx2 * cm, // G4double pDy2, G4double pDx3, G4double pDx4,
885  lg_Alp1 * rad // G4double pAlp2 //
886  );
887 
888  block_solid = new G4DisplacedSolid(G4String(GetName() + "_displaced"),
889  block_solid, 0, //
890  G4ThreeVector( //
891  tan(g_tower.pTheta * rad) * cos(g_tower.pPhi * rad), //
892  tan(g_tower.pTheta * rad) * sin(g_tower.pPhi * rad), //
893  1) * // G4ThreeVector
894  -(g_tower.pDz) *
895  cm //
896  + G4ThreeVector(shift_xcenter * cm, shift_ycenter * cm, 0) // shit in subtower direction
897  + G4ThreeVector(0, 0, -0.5 * g_tower.LightguideHeight * cm) //shift in the light guide height
898  );
899 
900  G4Material* cylinder_mat = G4Material::GetMaterial(
901  g_tower.LightguideMaterial);
902  assert(cylinder_mat);
903 
904  G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, cylinder_mat,
905  G4String(G4String(GetName()) + string("_Tower") + sTowerID), 0, 0,
906  step_limits);
907 
908  G4VisAttributes* VisAtt = new G4VisAttributes();
909  PHG4Utils::SetColour(VisAtt, g_tower.LightguideMaterial);
910  // VisAtt->SetColor(.3, .3, .3, .3);
911  VisAtt->SetVisibility(
912  _geom->is_azimuthal_seg_visible() or _geom->is_virualize_fiber());
913  VisAtt->SetForceSolid(not _geom->is_virualize_fiber());
914  block_logic->SetVisAttributes(VisAtt);
915 
916  return block_logic;
917 }