ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4SectorConstructor.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4SectorConstructor.cc
1 
11 
12 #include <g4main/PHG4DisplayAction.h> // for PHG4DisplayAction
13 #include <g4main/PHG4Subsystem.h> // for PHG4Subsystem
14 #include <g4main/PHG4Detector.h>
15 
16 #include <Geant4/G4Box.hh>
17 #include <Geant4/G4DisplacedSolid.hh> // for G4DisplacedSolid
18 #include <Geant4/G4Exception.hh> // for G4Exception
19 #include <Geant4/G4ExceptionSeverity.hh> // for FatalException, JustWarning
20 #include <Geant4/G4IntersectionSolid.hh>
21 #include <Geant4/G4LogicalVolume.hh>
22 #include <Geant4/G4Material.hh>
23 #include <Geant4/G4MaterialTable.hh> // for G4MaterialTable
24 #include <Geant4/G4PVPlacement.hh>
25 #include <Geant4/G4PhysicalConstants.hh> // for pi
26 #include <Geant4/G4Sphere.hh>
27 #include <Geant4/G4SystemOfUnits.hh> // for cm, um, perCent
28 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
29 #include <Geant4/G4Transform3D.hh> // for G4Transform3D, G4RotateX3D
30 #include <Geant4/G4Tubs.hh>
31 #include <Geant4/G4Types.hh> // for G4int
32 
33 #include <algorithm> // for max
34 #include <cassert>
35 #include <cmath>
36 #include <iostream>
37 #include <sstream>
38 #include <climits>
39 
40 using namespace PHG4Sector;
41 using namespace std;
42 
44  : overlapcheck_sector(false)
45  , name_base(name)
46  , m_DisplayAction(dynamic_cast<PHG4SectorDisplayAction *>(subsys->GetDisplayAction()))
47  , m_Verbosity(0)
48 {
49 }
50 
52 {
53  // geometry checks
54  if (geom.get_total_thickness() == 0)
55  {
57  (string("PHG4SectorConstructor::Construct_Sectors::") + (name_base)).c_str(),
58  __FILE__, FatalException,
59  " detector configured with zero thickness!");
60  }
61 
63  {
65  (string("PHG4SectorConstructor::Construct_Sectors::") + (name_base)).c_str(),
66  __FILE__, FatalException, "min_polar_angle = max_polar_angle!");
67  }
69  {
71  (string("PHG4SectorConstructor::Construct_Sectors::") + (name_base)).c_str(),
72  __FILE__, JustWarning,
73  "min and max polar angle got reversed. Correcting them.");
74 
75  const double t = geom.get_max_polar_angle();
78  }
80  {
82  (string("PHG4SectorConstructor::Construct_Sectors::") + (name_base)).c_str(),
83  __FILE__, FatalException,
84  "can NOT use flat edge for single or double sector detector!");
85  }
86 
87  const G4Transform3D transform_Det_to_Hall =
90 
91  const G4Transform3D transform_Hall_to_Det(transform_Det_to_Hall.inverse());
92 
93  // during GDML export, numerical value may change at the large digit and go beyond the 0 or pi limit.
94  // therefore recess 0/pi by a small amount to avoid such problem
95  static const double epsilon = numeric_limits<float>::epsilon();
96  const double sph_min_polar_angle =
98  const double sph_max_polar_angle =
100 
101  G4VSolid *SecConeBoundary_Hall = new G4Sphere("SecConeBoundary_Hall", //
102  geom.get_normal_start(), geom.get_max_R(), // G4double pRmin, G4double pRmax,
103  pi / 2 - pi / geom.get_N_Sector(), 2 * pi / geom.get_N_Sector(), // G4double pSPhi, G4double pDPhi,
104  sph_min_polar_angle, sph_max_polar_angle - sph_min_polar_angle //G4double pSTheta, G4double pDTheta
105  );
106 
107  G4VSolid *SecConeBoundary_Det = new G4DisplacedSolid("SecConeBoundary_Det",
108  SecConeBoundary_Hall, transform_Hall_to_Det);
109 
110  G4VSolid *Boundary_Det = SecConeBoundary_Det;
111 
113  {
114  // build a flat edge
115 
116  const double sph_min_polar_size =
118  const double sph_max_polar_size =
120 
121  G4VSolid *BoxBoundary_Det = new G4Box("BoxBoundary_Det",
122  geom.get_max_R(), (sph_min_polar_size + sph_max_polar_size) / 2,
124  G4VSolid *BoxBoundary_Det_Place = new G4DisplacedSolid(
125  "BoxBoundary_Det_Place", BoxBoundary_Det, 0,
126  G4ThreeVector(0, (sph_max_polar_size - sph_min_polar_size) / 2, 0));
127 
128  Boundary_Det = new G4IntersectionSolid("Boundary_Det",
129  BoxBoundary_Det_Place, SecConeBoundary_Det);
130  }
131 
132  G4VSolid *DetectorBox_Det = Construct_Sectors_Plane("DetectorBox_Det",
134  Boundary_Det);
135 
137 
138  G4LogicalVolume *DetectorLog_Det = new G4LogicalVolume(DetectorBox_Det, //
139  p_mat, name_base + "_Log");
140  RegisterLogicalVolume(DetectorLog_Det);
141 
142  for (G4int sec = 0; sec < geom.get_N_Sector(); sec++)
143  {
145  new G4PVPlacement(
146  G4RotateZ3D(2 * pi / geom.get_N_Sector() * sec) * transform_Det_to_Hall, DetectorLog_Det,
147  name_base + "_Physical", WorldLog, false, sec, overlapcheck_sector));
148  }
149 
150  // construct layers
151  double z_start = -geom.get_total_thickness() / 2;
152 
153  for (Sector_Geometry::t_layer_list::const_iterator it =
154  geom.layer_list.begin();
155  it != geom.layer_list.end(); ++it)
156  {
157  const Layer &l = (*it);
158 
159  if (l.percentage_filled > 100. or l.percentage_filled < 0)
160  {
161  stringstream s;
162  s << name_base << " have invalid layer ";
163  s << l.name << " with percentage_filled =" << l.percentage_filled;
164 
165  G4Exception(
166  (string("PHG4SectorConstructor::Construct_Sectors::") + (name_base)).c_str(), __FILE__, FatalException,
167  s.str().c_str());
168  }
169 
170  const string layer_name = name_base + "_" + l.name;
171 
172  G4VSolid *LayerSol_Det = Construct_Sectors_Plane(layer_name + "_Sol",
173  z_start, l.depth * l.percentage_filled * perCent, Boundary_Det);
174 
175  G4LogicalVolume *LayerLog_Det = new G4LogicalVolume(LayerSol_Det, //
176  PHG4Detector::GetDetectorMaterial(l.material), layer_name + "_Log");
177  RegisterLogicalVolume(LayerLog_Det);
178 
180  new G4PVPlacement(0, G4ThreeVector(), LayerLog_Det,
181  layer_name + "_Physical", DetectorLog_Det, false, 0, overlapcheck_sector),
182  l.active);
183 
184  z_start += l.depth;
185  }
186 
187  if (abs(z_start - geom.get_total_thickness() / 2) > 1 * um)
188  {
189  stringstream s;
190  s << name_base << " - accumulated thickness = "
191  << (z_start + geom.get_total_thickness() / 2) / um
192  << " um expected thickness = " << geom.get_total_thickness() / um
193  << " um";
194  G4Exception(
195  (string("PHG4SectorConstructor::Construct_Sectors::") + (name_base)).c_str(),
196  __FILE__, FatalException, s.str().c_str());
197  }
198 
199  m_DisplayAction->AddVolume(DetectorLog_Det, "DetectorBox");
200  if (Verbosity() > 1)
201  {
202  std::cout << "PHG4SectorConstructor::Construct_Sectors::" << name_base
203  << " - total thickness = " << geom.get_total_thickness() / cm << " cm"
204  << std::endl;
205  std::cout << "PHG4SectorConstructor::Construct_Sectors::" << name_base << " - "
206  << map_log_vol.size() << " logical volume constructed" << std::endl;
207  std::cout << "PHG4SectorConstructor::Construct_Sectors::" << name_base << " - "
208  << map_phy_vol.size() << " physical volume constructed; "
209  << map_active_phy_vol.size() << " is active." << std::endl;
210  }
211 }
212 
213 G4VSolid *
215  const std::string &name, //
216  const double start_z, //
217  const double thickness, //
218  G4VSolid *SecConeBoundary_Det //
219 )
220 {
221  assert(SecConeBoundary_Det);
222 
223  G4VSolid *Sol_Raw = new G4Tubs(name + "_Raw", //const G4String& pName,
224  0, // G4double pRMin,
225  geom.get_max_R(), // G4double pRMax,
226  thickness / 2, // G4double pDz,
227  0, // G4double pSPhi,
228  2 * pi // G4double pDPhi
229  );
230 
231  G4VSolid *Sol_Place = new G4DisplacedSolid(name + "_Place", Sol_Raw, 0,
232  G4ThreeVector(0, 0, start_z + thickness / 2));
233 
234  G4VSolid *Sol = new G4IntersectionSolid(name.c_str(), Sol_Place,
235  SecConeBoundary_Det);
236 
237  return Sol;
238  // return Sol_Place;
239 }
240 
242 {
243  N_Sector = 8;
244  material = "G4_AIR";
245 
247  normal_polar_angle = 0;
248 
250  min_polar_angle = .1;
251 
253  max_polar_angle = .3;
254 
256  normal_start = 305 * cm;
257 
259 
261 }
262 
265 {
266  if (!v)
267  {
268  std::cout
269  << "PHG4SectorConstructor::RegisterVolume - Error - invalid volume!"
270  << std::endl;
271  return v;
272  }
273  if (map_log_vol.find(v->GetName()) != map_log_vol.end())
274  {
275  std::cout << "PHG4SectorConstructor::RegisterVolume - Warning - replacing "
276  << v->GetName() << std::endl;
277  }
278 
279  map_log_vol[v->GetName()] = v;
280 
281  return v;
282 }
283 
286  const bool active)
287 {
288  if (!v)
289  {
290  std::cout
291  << "PHG4SectorConstructor::RegisterPhysicalVolume - Error - invalid volume!"
292  << std::endl;
293  return v;
294  }
295 
296  phy_vol_idx_t id(v->GetName(), v->GetCopyNo());
297 
298  if (map_phy_vol.find(id) != map_phy_vol.end())
299  {
300  std::cout
301  << "PHG4SectorConstructor::RegisterPhysicalVolume - Warning - replacing "
302  << v->GetName() << "[" << v->GetCopyNo() << "]" << std::endl;
303  }
304 
305  map_phy_vol[id] = v;
306 
307  if (active)
308  map_active_phy_vol[id] = v;
309 
310  return v;
311 }
312 
313 double
315 {
316  double sum = 0;
317  for (t_layer_list::const_iterator it = layer_list.begin();
318  it != layer_list.end(); ++it)
319  {
320  sum += (*it).depth;
321  }
322  return sum;
323 }
324 
325 double
327 {
328  // Geometry check
329  assert(abs(min_polar_angle - normal_polar_angle) < pi / 2);
330  assert(abs(max_polar_angle - normal_polar_angle) < pi / 2);
331 
332  if (N_Sector <= 2)
333  {
334  // Geometry check
335  if (cos(min_polar_angle + normal_polar_angle) <= 0)
336  {
337  stringstream s;
338  s << "Geometry check failed. " << endl;
339  s << "normal_polar_angle = " << normal_polar_angle << endl;
340  s << "min_polar_angle = " << min_polar_angle << endl;
341  s << "max_polar_angle = " << max_polar_angle << endl;
342  s << "cos(min_polar_angle + normal_polar_angle) = "
343  << cos(min_polar_angle + normal_polar_angle) << endl;
344 
345  G4Exception("Sector_Geometry::get_max_R", __FILE__, FatalException,
346  s.str().c_str());
347  }
348  if (cos(max_polar_angle + normal_polar_angle) <= 0)
349  {
350  stringstream s;
351  s << "Geometry check failed. " << endl;
352  s << "normal_polar_angle = " << normal_polar_angle << endl;
353  s << "min_polar_angle = " << min_polar_angle << endl;
354  s << "max_polar_angle = " << max_polar_angle << endl;
355  s << "cos(max_polar_angle + normal_polar_angle) = "
356  << cos(max_polar_angle + normal_polar_angle) << endl;
357 
358  G4Exception("Sector_Geometry::get_max_R", __FILE__, FatalException,
359  s.str().c_str());
360  }
361 
362  const double max_tan_angle = std::max(
365 
366  return (get_normal_start() + get_total_thickness()) * sqrt(1 + max_tan_angle * max_tan_angle);
367  }
368  else
369  {
370  const double max_angle = std::max(
373 
374  return (get_normal_start() + get_total_thickness()) * sqrt(1 + pow(tan(max_angle), 2) + pow(tan(2 * pi / N_Sector), 2)) * 2;
375  }
376 }
377 
381 void Sector_Geometry::AddLayers_DriftVol_COMPASS(const double drift_vol_thickness)
382 {
383  // (drift chamber) is enclosed by two Mylarr [52] cathode foils of 25um thickness,
384  // coated with about 10um of graphite,
385 
386  AddLayer("EntranceWindow", "G4_MYLAR", 25 * um, false, 100);
387  AddLayer("Cathode", "G4_GRAPHITE", 10 * um, false, 100);
388  AddLayer("DrftVol", material, drift_vol_thickness, true, 100);
389 }
390 
395 void Sector_Geometry::AddLayers_HBD_GEM(const int n_GEM_layers)
396 {
397  // Internal HBD structure
398  // From doi:10.1016/j.nima.2011.04.015
399  // Component Material X0 (cm) Thickness (cm) Area (%) Rad. Length (%)
400  // Mesh SS 1.67 0.003 11.5 0.021 <- not used for GEMs trackers
401  // AddLayer("Mesh", "Steel",
402  // 0.003 * cm, false, 11.5);
403 
404  // // GEM frames FR4 17.1 0.15x4 6.5 0.228 <- not used for GEMs trackers
405  // AddLayer("Frame0", "G10",
406  // 0.15 * cm, false, 6.5);
407 
408  for (int gem = 1; gem <= n_GEM_layers; gem++)
409  {
410  stringstream sid;
411  sid << gem;
412 
413  // GEM Copper 1.43 0.0005x6 64 0.134
414  AddLayer(G4String("GEMFrontCu") + G4String(sid.str()), "G4_Cu",
415  0.0005 * cm, false, 64);
416 
417  // GEM Kapton 28.6 0.005x3 64 0.034
418  AddLayer(G4String("GEMKapton") + G4String(sid.str()), "G4_KAPTON",
419  0.005 * cm, false, 64);
420 
421  // GEM Copper 1.43 0.0005x6 64 0.134
422  AddLayer(G4String("GEMBackCu") + G4String(sid.str()), "G4_Cu",
423  0.0005 * cm, false, 64);
424 
425  // GEM frames FR4 17.1 0.15x4 6.5 0.228
426  AddLayer(G4String("Frame") + G4String(sid.str()), "G10", 0.15 * cm, false,
427  6.5);
428  }
429 
430  // PCB Kapton 28.6 0.005 100 0.017
431  AddLayer(G4String("PCBKapton"), "G4_KAPTON", 0.005 * cm, false, 100);
432 
433  // PCB Copper 1.43 0.0005 80 0.028
434  AddLayer(G4String("PCBCu"), "G4_Cu", 0.0005 * cm, false, 80);
435 
436  // Facesheet FR4 17.1 0.025x2 100 0.292
437  AddLayer("Facesheet", "G10", 0.025 * 2 * cm, false, 100);
438 
439  // Panel core Honeycomb 8170 1.905 100 0.023 <- very thin-X0 stuff, ignore
440 
441  // Total vessel 0.82
442  // Readout
443 }
444 
450 {
451  // Readout
452  // Readout board FR4/copper 17.1/1.43 0.05/0.001 100 0.367
453  AddLayer(G4String("ReadoutFR4"), "G10", 0.05 * cm, false, 100);
454  AddLayer(G4String("ReadoutCu"), "G4_Cu", 0.001 * cm, false, 100);
455 
456  // Preamps + sockets Copper 1.43 0.0005 100 0.66
457  // Total readout 1.03
458  AddLayer(G4String("SocketsCu"), "G4_Cu", 0.0005 * cm, false, 100);
459 }
460 
465 void Sector_Geometry::AddLayers_AeroGel_ePHENIX(const double radiator_length,
466  const double expansion_length, std::string radiator)
467 {
468  if (radiator == "Default")
469  {
470  radiator = "ePHENIX_AeroGel";
471  }
472 
473  // Readout board FR4/copper 17.1/1.43 0.05/0.001 100 0.367
474  AddLayer("EntranceWindow", "G10", 0.05 * cm, false, 100);
475  AddLayer("AeroGel", radiator, radiator_length, false, 100);
476  AddLayer("ExpansionVol", "G4_AIR", expansion_length, false, 100);
477 
478  //Some readout
479  AddLayer(G4String("ReadoutFR4"), "G10", 0.05 * cm, false, 100);
480  AddLayer(G4String("ReadoutCu"), "G4_Cu", 0.001 * cm, false, 100);
481  AddLayer(G4String("SocketsCu"), "G4_Cu", 0.0005 * cm, false, 100);
482 }