EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4BarrelEcalDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4BarrelEcalDetector.cc
1 #include "PHG4BarrelEcalDetector.h"
2 #include "PHG4BarrelEcalDisplayAction.h"
3 
4 #include <phparameter/PHParameters.h>
5 
6 #include <TSystem.h>
7 #include <g4main/PHG4Detector.h> // for PHG4Detector
8 #include <g4main/PHG4DisplayAction.h> // for PHG4DisplayAction
9 #include <g4main/PHG4Subsystem.h>
10 #include <phool/recoConsts.h>
11 #include <Geant4/G4Box.hh>
12 #include <Geant4/G4Cons.hh>
13 #include <Geant4/G4CutTubs.hh>
14 #include <Geant4/G4DisplacedSolid.hh>
15 #include <Geant4/G4LogicalVolume.hh>
16 #include <Geant4/G4Material.hh>
17 #include <Geant4/G4MaterialPropertiesTable.hh> // for G4MaterialProperties...
18 #include <Geant4/G4MaterialPropertyVector.hh> // for G4MaterialPropertyVector
19 #include <Geant4/G4PVPlacement.hh>
20 #include <Geant4/G4RotationMatrix.hh> // for G4RotationMatrix
21 #include <Geant4/G4SubtractionSolid.hh>
22 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
23 #include <Geant4/G4Transform3D.hh> // for G4Transform3D
24 #include <Geant4/G4Tubs.hh>
25 #include <Geant4/G4Types.hh> // for G4double, G4int
26 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
27 #include <Geant4/G4NistManager.hh>
28 
29 #include <g4gdml/PHG4GDMLConfig.hh>
30 #include <g4gdml/PHG4GDMLUtility.hh>
31 
32 #include <cmath>
33 #include <cstdlib>
34 #include <fstream>
35 #include <iostream>
36 #include <sstream>
37 #include <utility> // for pair, make_pair
38 
39 using namespace std;
40 
41 class G4VSolid;
42 class PHCompositeNode;
43 
44 //_______________________________________________________________________
45 PHG4BarrelEcalDetector::PHG4BarrelEcalDetector(PHG4Subsystem* subsys, PHCompositeNode* Node, PHParameters* parameters, const std::string& dnam)
46  : PHG4Detector(subsys, Node, dnam)
47  , m_DisplayAction(dynamic_cast<PHG4BarrelEcalDisplayAction*>(subsys->GetDisplayAction()))
48  , m_Params(parameters)
49  , m_ActiveFlag(m_Params->get_int_param("active"))
50  , m_AbsorberActiveFlag(m_Params->get_int_param("absorberactive"))
51  , m_SupportActiveFlag(m_Params->get_int_param("supportactive"))
52  , m_TowerLogicNamePrefix("bcalTower")
53  , m_SuperDetector("NONE")
54 {
55  gdml_config = PHG4GDMLUtility::GetOrMakeConfigNode(Node);
56  assert(gdml_config);
57 }
58 //_______________________________________________________________________
59 int PHG4BarrelEcalDetector::IsInBarrelEcal(G4VPhysicalVolume* volume) const
60 {
61  G4LogicalVolume* mylogvol = volume->GetLogicalVolume();
62  if (m_ActiveFlag)
63  {
64  if (m_ScintiLogicalVolSet.find(mylogvol) != m_ScintiLogicalVolSet.end())
65  {
66  return 1;
67  }
68  }
69 
70  if (m_AbsorberActiveFlag)
71  {
72  if (m_AbsorberLogicalVolSet.find(mylogvol) != m_AbsorberLogicalVolSet.end())
73  {
74  return -1;
75  }
76  }
77 
78  if (m_SupportActiveFlag)
79  {
80  if (m_SupportLogicalVolSet.find(mylogvol) != m_SupportLogicalVolSet.end())
81  {
82  return -2;
83  }
84  }
85  return 0;
86 }
87 
88 //_______________________________________________________________________
89 void PHG4BarrelEcalDetector::ConstructMe(G4LogicalVolume* logicWorld)
90 {
91  if (Verbosity() > 0)
92  {
93  std::cout << "PHG4BarrelEcalDetector: Begin Construction" << std::endl;
94  }
95 
96  if (m_Params->get_string_param("mapping_file").empty())
97  {
98  std::cout << "ERROR in PHG4BarrelEcalDetector: No mapping file specified. Abort detector construction." << std::endl;
99  std::cout << "Please run set_string_param(\"mapping_file\", std::string filename ) first." << std::endl;
100  gSystem->Exit(1);
101  }
102 
103  ParseParametersFromTable();
104 
105  G4double Radius = m_Params->get_double_param("radius") * cm;
106  G4double tower_length = m_Params->get_double_param("tower_length") * cm;
107  G4double becal_length = m_Params->get_double_param("becal_length") * cm;
108  G4double CenterZ_Shift = m_Params->get_double_param("CenterZ_Shift") * cm;
109  G4double cone1_h = m_Params->get_double_param("cone1_h") * cm;
110  G4double cone1_dz = m_Params->get_double_param("cone1_dz") * cm;
111  G4double cone2_h = m_Params->get_double_param("cone2_h") * cm;
112  G4double cone2_dz = m_Params->get_double_param("cone2_dz") * cm;
113 
114  silicon_width_half = m_Params->get_double_param("silicon_width_half") * cm;
115  kapton_width_half = m_Params->get_double_param("kapton_width_half") * cm;
116  SIO2_width_half = m_Params->get_double_param("SIO2_width_half") * cm;
117  Carbon_width_half = m_Params->get_double_param("Carbon_width_half") * cm;
118  support_length = m_Params->get_double_param("support_length") * cm;
119 
120  G4double max_radius = Radius + tower_length + 2 * silicon_width_half + 2 * kapton_width_half + 2 * SIO2_width_half + 2 * Carbon_width_half + support_length + 8 * overlap;
121 
122  //std::cout << Radius << " " << max_radius << "====================================" << std::endl;
123 
124  G4double pos_x1 = 0 * cm;
125  G4double pos_y1 = 0 * cm;
126  G4double pos_z1 = 0 * cm;
127 
128  G4Tubs* cylinder_solid1 = new G4Tubs("BCAL_SOLID1",
129  Radius, max_radius,
130  becal_length / 2, 0, 2 * M_PI);
131 
132  G4Tubs* cylinder_solid2 = new G4Tubs("BCAL_SOLID2",
133  Radius - 1, max_radius + 1,
134  abs(CenterZ_Shift), 0, 2 * M_PI);
135 
136  G4ThreeVector shift_cs2 = G4ThreeVector(0, 0, becal_length / 2);
137 
138  G4VSolid* cylinder_solid3 = new G4SubtractionSolid("BCAL_SOLID3", cylinder_solid1, cylinder_solid2, 0, shift_cs2);
139 
140  G4Cons* cone1 = new G4Cons("cone1",
141  Radius - 1, Radius - 1,
142  Radius - 1, Radius + cone1_h,
143  cone1_dz, 0, 2 * M_PI);
144 
145  G4ThreeVector shift_cone1 = G4ThreeVector(0, 0, becal_length / 2 - abs(CenterZ_Shift) - cone1_dz);
146 
147  G4VSolid* cylinder_solid4 = new G4SubtractionSolid("BCAL_SOLID4", cylinder_solid3, cone1, 0, shift_cone1);
148 
149  G4Cons* cone2 = new G4Cons("cone2",
150  Radius - 1, Radius + cone2_h,
151  Radius - 1, Radius - 1,
152  cone2_dz, 0, 2 * M_PI);
153 
154  G4ThreeVector shift_cone2 = G4ThreeVector(0, 0, -becal_length / 2 + cone2_dz);
155 
156  G4VSolid* cylinder_solid = new G4SubtractionSolid("BCAL_SOLID", cylinder_solid4, cone2, 0, shift_cone2);
157 
158  G4Material* cylinder_mat = GetDetectorMaterial("G4_AIR");
159  assert(cylinder_mat);
160 
161  G4LogicalVolume* cylinder_logic = new G4LogicalVolume(cylinder_solid, cylinder_mat,
162  "BCAL_SOLID", 0, 0, 0);
163 
164  m_DisplayAction->AddVolume(cylinder_logic, "BCalCylinder");
165 
166  std::string name_envelope = m_TowerLogicNamePrefix + "_envelope";
167 
168  G4PVPlacement* phys_envelope =
169  new G4PVPlacement(0, G4ThreeVector(pos_x1, pos_y1, pos_z1), cylinder_logic, name_envelope,
170  logicWorld, false, 0, OverlapCheck());
171 
172  gdml_config->exclude_physical_vol(phys_envelope);
173  PlaceTower(cylinder_logic);
174 
175  return;
176 }
177 
178 int PHG4BarrelEcalDetector::PlaceTower(G4LogicalVolume* sec)
179 {
180  /* Loop over all tower positions in vector and place tower */
181  for (std::map<std::string, towerposition>::iterator iterator = m_TowerPostionMap.begin(); iterator != m_TowerPostionMap.end(); ++iterator)
182  {
183  if (Verbosity() > 0)
184  {
185  std::cout << "PHG4BarrelEcalDetector: Place tower " << iterator->first
186  << " idx_j = " << iterator->second.idx_j << ", idx_k = " << iterator->second.idx_k << std::endl;
187  }
188 
189  int copyno = (iterator->second.idx_j << 16) + iterator->second.idx_k;
190 
191  G4LogicalVolume* block_logic = ConstructTower(iterator);
192  m_DisplayAction->AddVolume(block_logic, iterator->first);
193  if (iterator->second.idx_k % 2 == 0)
194  {
195  m_DisplayAction->AddVolume(block_logic, "Block1");
196  }
197  else
198  {
199  m_DisplayAction->AddVolume(block_logic, "Block2");
200  }
201 
202  G4LogicalVolume* glass_logic = ConstructGlass(iterator);
203  if (iterator->second.idx_k % 2 == 0)
204  {
205  m_DisplayAction->AddVolume(glass_logic, "Block1");
206  }
207  else
208  {
209  m_DisplayAction->AddVolume(glass_logic, "Block2");
210  }
211 
212  G4RotationMatrix becal_rotm;
213  becal_rotm.rotateY(iterator->second.rotx);
214  becal_rotm.rotateY(iterator->second.roty);
215  becal_rotm.rotateZ(iterator->second.rotz);
216 
217  new G4PVPlacement(G4Transform3D(becal_rotm, G4ThreeVector(iterator->second.centerx, iterator->second.centery, iterator->second.centerz)),
218  block_logic,
219  G4String(string(iterator->first) + string("_TT")),
220  sec,
221  0, copyno, OverlapCheck());
222 
223  G4double posx_glass = iterator->second.centerx + th / 2 * abs(cos(iterator->second.roty - M_PI_2)) * cos(iterator->second.rotz);
224  G4double posy_glass = iterator->second.centery + th / 2 * abs(cos(iterator->second.roty - M_PI_2)) * sin(iterator->second.rotz);
225  G4double posz_glass = iterator->second.centerz + th * sin(iterator->second.roty - M_PI_2);
226 
227  new G4PVPlacement(G4Transform3D(becal_rotm, G4ThreeVector(posx_glass, posy_glass, posz_glass)),
228  glass_logic,
229  G4String(string(iterator->first) + string("_Glass")),
230  sec,
231  0, copyno, OverlapCheck());
232 
233  //=================Silicon
234  G4LogicalVolume* block_silicon = ConstructSi(iterator);
235  m_DisplayAction->AddVolume(block_silicon, "Si");
236 
237  G4double pTheta = iterator->second.pTheta;
238  G4double theta = iterator->second.roty - M_PI_2;
239  G4double len = (tower_length / 2) + silicon_width_half + overlap;
240  G4double sci_sr = len * tan(pTheta);
241  G4double sci_sz = len;
242  G4double sci_mag = sqrt(sci_sr * sci_sr + sci_sz * sci_sz);
243 
244  G4double posz_si = iterator->second.centerz - sci_mag * sin(theta + pTheta);
245  G4double posy_si = iterator->second.centery + sci_mag * cos(theta + pTheta) * sin(iterator->second.rotz);
246  G4double posx_si = iterator->second.centerx + sci_mag * cos(theta + pTheta) * cos(iterator->second.rotz);
247  new G4PVPlacement(G4Transform3D(becal_rotm, G4ThreeVector(posx_si, posy_si, posz_si)),
248  block_silicon,
249  G4String(string(iterator->first) + string("_Si")),
250  sec,
251  0, copyno, OverlapCheck());
252 
253  //=================Kapton
254  G4LogicalVolume* block_kapton = ConstructKapton(iterator);
255  m_DisplayAction->AddVolume(block_kapton, "Kapton");
256 
257  len += kapton_width_half + overlap + silicon_width_half;
258 
259  G4double kapton_sr = len * tan(pTheta);
260  G4double kapton_sz = len;
261  G4double kapton_mag = sqrt(kapton_sr * kapton_sr + kapton_sz * kapton_sz);
262 
263  G4double posz_kapton = iterator->second.centerz - kapton_mag * sin(theta + pTheta);
264  G4double posy_kapton = iterator->second.centery + kapton_mag * cos(theta + pTheta) * sin(iterator->second.rotz);
265  G4double posx_kapton = iterator->second.centerx + kapton_mag * cos(theta + pTheta) * cos(iterator->second.rotz);
266  new G4PVPlacement(G4Transform3D(becal_rotm, G4ThreeVector(posx_kapton, posy_kapton, posz_kapton)),
267  block_kapton,
268  G4String(string(iterator->first) + string("_Kapton")),
269  sec,
270  0, copyno, OverlapCheck());
271 
272  //=================SIO2
273  G4LogicalVolume* block_SIO2 = ConstructSIO2(iterator);
274  m_DisplayAction->AddVolume(block_SIO2, "SIO2");
275 
276  len += SIO2_width_half + overlap + kapton_width_half;
277 
278  G4double SIO2_sr = len * tan(pTheta);
279  G4double SIO2_sz = len;
280  G4double SIO2_mag = sqrt(SIO2_sr * SIO2_sr + SIO2_sz * SIO2_sz);
281 
282  G4double posz_SIO2 = iterator->second.centerz - SIO2_mag * sin(theta + pTheta);
283  G4double posy_SIO2 = iterator->second.centery + SIO2_mag * cos(theta + pTheta) * sin(iterator->second.rotz);
284  G4double posx_SIO2 = iterator->second.centerx + SIO2_mag * cos(theta + pTheta) * cos(iterator->second.rotz);
285  new G4PVPlacement(G4Transform3D(becal_rotm, G4ThreeVector(posx_SIO2, posy_SIO2, posz_SIO2)),
286  block_SIO2,
287  G4String(string(iterator->first) + string("SIO2")),
288  sec,
289  0, copyno, OverlapCheck());
290 
291  //=================Carbon
292  G4LogicalVolume* block_Carbon = ConstructCarbon(iterator);
293 
294  if (iterator->second.idx_k % 2 == 0)
295  {
296  m_DisplayAction->AddVolume(block_Carbon, "Block1");
297  }
298  else
299  {
300  m_DisplayAction->AddVolume(block_Carbon, "Block2");
301  }
302 
303  len += Carbon_width_half + overlap + SIO2_width_half;
304 
305  G4double Carbon_sr = len * tan(pTheta);
306  G4double Carbon_sz = len;
307  G4double Carbon_mag = sqrt(Carbon_sr * Carbon_sr + Carbon_sz * Carbon_sz);
308 
309  G4double posz_Carbon = iterator->second.centerz - Carbon_mag * sin(theta + pTheta);
310  G4double posy_Carbon = iterator->second.centery + Carbon_mag * cos(theta + pTheta) * sin(iterator->second.rotz);
311  G4double posx_Carbon = iterator->second.centerx + Carbon_mag * cos(theta + pTheta) * cos(iterator->second.rotz);
312  new G4PVPlacement(G4Transform3D(becal_rotm, G4ThreeVector(posx_Carbon, posy_Carbon, posz_Carbon)),
313  block_Carbon,
314  G4String(string(iterator->first) + string("Carbon")),
315  sec,
316  0, copyno, OverlapCheck());
317  }
318  return 0;
319 }
320 
321 int PHG4BarrelEcalDetector::ParseParametersFromTable()
322 {
323  /* Open the datafile, if it won't open return an error */
324  std::ifstream istream_mapping;
325  istream_mapping.open(m_Params->get_string_param("mapping_file"));
326  if (!istream_mapping.is_open())
327  {
328  std::cout << "ERROR in PHG4BarrelEcalDetector: Failed to open mapping file " << m_Params->get_string_param("mapping_file") << std::endl;
329  gSystem->Exit(1);
330  }
331 
332  /* loop over lines in file */
333  std::string line_mapping;
334  while (getline(istream_mapping, line_mapping))
335  {
336  std::istringstream iss(line_mapping);
337 
338  if (line_mapping.find("BECALtower ") != string::npos)
339  {
340  unsigned idphi_j, ideta_k;
341  G4double cx, cy, cz;
342  G4double rot_z, rot_y, rot_x;
343  G4double size_z, size_y1, size_x1, size_y2, size_x2, p_Theta;
344  std::string dummys;
345 
346  if (!(iss >> dummys >> ideta_k >> idphi_j >> size_x1 >> size_x2 >> size_y1 >> size_y2 >> size_z >> p_Theta >> cx >> cy >> cz >> rot_x >> rot_y >> rot_z))
347  {
348  std::cout << "ERROR in PHG4BarrelEcalDetector: Failed to read line in mapping file " << m_Params->get_string_param("mapping_file") << std::endl;
349  gSystem->Exit(1);
350  }
351 
352  /* Construct unique name for tower */
353  /* Mapping file uses cm, this class uses mm for length */
354  std::ostringstream towername;
355  towername.str("");
356  towername << m_TowerLogicNamePrefix << "_j_" << idphi_j << "_k_" << ideta_k;
357 
358  /* insert tower into tower map */
359  towerposition tower_new;
360  tower_new.sizex1 = size_x1 * cm;
361  tower_new.sizex2 = size_x2 * cm;
362  tower_new.sizey1 = size_y1 * cm;
363  tower_new.sizey2 = size_y2 * cm;
364  tower_new.sizez = size_z * cm;
365  tower_new.pTheta = p_Theta;
366  tower_new.centerx = cx * cm;
367  tower_new.centery = cy * cm;
368  tower_new.centerz = cz * cm;
369  tower_new.rotx = rot_x;
370  tower_new.roty = rot_y;
371  tower_new.rotz = rot_z;
372  tower_new.idx_j = idphi_j;
373  tower_new.idx_k = ideta_k;
374  m_TowerPostionMap.insert(make_pair(towername.str(), tower_new));
375  }
376  else
377  {
378  /* If this line is not a comment and not a tower, save parameter as string / value. */
379  string parname;
380  G4double parval;
381 
382  /* read string- break if error */
383  if (!(iss >> parname >> parval))
384  {
385  cout << "ERROR in PHG4CrystalCalorimeterDetector: Failed to read line in mapping file " << m_Params->get_string_param("mappingtower") << endl;
386  gSystem->Exit(1);
387  }
388 
389  m_GlobalParameterMap.insert(make_pair(parname, parval));
390 
391  /* Update member variables for global parameters based on parsed parameter file */
392 
393  std::map<string, G4double>::iterator parit;
394 
395  parit = m_GlobalParameterMap.find("CenterZ_Shift");
396  if (parit != m_GlobalParameterMap.end())
397  {
398  m_Params->set_double_param("CenterZ_Shift", parit->second); // in cm
399  }
400  parit = m_GlobalParameterMap.find("radius");
401  if (parit != m_GlobalParameterMap.end())
402  {
403  m_Params->set_double_param("radius", parit->second); // in cm
404  }
405  parit = m_GlobalParameterMap.find("tower_length");
406  if (parit != m_GlobalParameterMap.end())
407  {
408  m_Params->set_double_param("tower_length", parit->second); // in cm
409  }
410  parit = m_GlobalParameterMap.find("becal_length");
411  if (parit != m_GlobalParameterMap.end())
412  {
413  m_Params->set_double_param("becal_length", parit->second); // in cm
414  }
415  parit = m_GlobalParameterMap.find("cone1_h");
416  if (parit != m_GlobalParameterMap.end())
417  {
418  m_Params->set_double_param("cone1_h", parit->second); // in cm
419  }
420  parit = m_GlobalParameterMap.find("cone1_dz");
421  if (parit != m_GlobalParameterMap.end())
422  {
423  m_Params->set_double_param("cone1_dz", parit->second); // in cm
424  }
425  parit = m_GlobalParameterMap.find("cone2_h");
426  if (parit != m_GlobalParameterMap.end())
427  {
428  m_Params->set_double_param("cone2_dz", parit->second); // in cm
429  }
430  if (parit != m_GlobalParameterMap.end())
431  {
432  m_Params->set_double_param("cone2_dz", parit->second); // in cm
433  }
434  parit = m_GlobalParameterMap.find("thickness_wall");
435  if (parit != m_GlobalParameterMap.end())
436  {
437  m_Params->set_double_param("thickness_wall", parit->second); // in cm
438  }
439  parit = m_GlobalParameterMap.find("silicon_width_half");
440  if (parit != m_GlobalParameterMap.end())
441  {
442  m_Params->set_double_param("silicon_width_half", parit->second); // in cm
443  }
444  parit = m_GlobalParameterMap.find("kapton_width_half");
445  if (parit != m_GlobalParameterMap.end())
446  {
447  m_Params->set_double_param("kapton_width_half", parit->second); // in cm
448  }
449  parit = m_GlobalParameterMap.find("SIO2_width_half");
450  if (parit != m_GlobalParameterMap.end())
451  {
452  m_Params->set_double_param("SIO2_width_half", parit->second); // in cm
453  }
454  parit = m_GlobalParameterMap.find("Carbon_width_half");
455  if (parit != m_GlobalParameterMap.end())
456  {
457  m_Params->set_double_param("Carbon_width_half", parit->second); // in cm
458  }
459  parit = m_GlobalParameterMap.find("support_length");
460  if (parit != m_GlobalParameterMap.end())
461  {
462  m_Params->set_double_param("support_length", parit->second); // in cm
463  }
464  }
465  }
466 
467  return 0;
468 }
469 
470 G4LogicalVolume*
471 PHG4BarrelEcalDetector::ConstructTower(std::map<std::string, towerposition>::iterator iterator)
472 {
473  G4Trap* block_tower = GetTowerTrap(iterator);
474  G4Trap* block_glass = GetGlassTrapSubtract(iterator);
475  G4ThreeVector shift = G4ThreeVector(-th * sin(iterator->second.roty - M_PI_2), 0, th / 2 * abs(cos(iterator->second.roty - M_PI_2)));
476  G4VSolid* block_solid = new G4SubtractionSolid(G4String(string(iterator->first) + string("_Envelope")), block_tower, block_glass, 0, shift);
477  G4Material* material_shell = GetCarbonFiber();
478  assert(material_shell);
479 
480  G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, material_shell,
481  G4String(string(iterator->first) + string("_Tower")), 0, 0,
482  nullptr);
483  m_AbsorberLogicalVolSet.insert(block_logic);
484  return block_logic;
485 }
486 
487 G4LogicalVolume*
488 PHG4BarrelEcalDetector::ConstructGlass(std::map<std::string, towerposition>::iterator iterator)
489 {
490  G4Trap* block_solid = GetGlassTrap(iterator);
491  G4Material* material_glass = GetSciGlass();
492  assert(material_glass);
493  G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, material_glass,
494  G4String(string(iterator->first) + string("_Glass")), 0, 0,
495  nullptr);
496  m_ScintiLogicalVolSet.insert(block_logic);
497  return block_logic;
498 }
499 
500 G4Material* PHG4BarrelEcalDetector::GetCarbonFiber()
501 {
502  static string matname = "CrystalCarbonFiber";
503  G4Material* carbonfiber = GetDetectorMaterial(matname, false); // false suppresses warning that material does not exist
504  if (!carbonfiber)
505  {
506  G4double density_carbon_fiber = 1.44 * g / cm3;
507  carbonfiber = new G4Material(matname, density_carbon_fiber, 1);
508  carbonfiber->AddElement(GetDetectorElement("C"), 1);
509  }
510  return carbonfiber;
511 }
512 
513 G4Material* PHG4BarrelEcalDetector::GetSciGlass()
514 {
515  static string matname = "sciglass";
516  G4Material* sciglass = GetDetectorMaterial(matname, false); // false suppresses warning that material does not exist
517  if (!sciglass)
518  {
519  G4double density;
520  G4int ncomponents;
521 
522  sciglass = new G4Material(matname, density = 4.22 * g / cm3, ncomponents = 4, kStateSolid);
523  sciglass->AddElement(GetDetectorElement("Ba"), 0.3875);
524  sciglass->AddElement(GetDetectorElement("Gd"), 0.2146);
525  sciglass->AddElement(GetDetectorElement("Si"), 0.1369);
526  sciglass->AddElement(GetDetectorElement("O"), 0.2610);
527  }
528 
529  return sciglass;
530 }
531 
532 G4Trap* PHG4BarrelEcalDetector::GetGlassTrap(std::map<std::string, towerposition>::iterator iterator)
533 {
534  G4double th = m_Params->get_double_param("thickness_wall") * cm;
535  G4double size_x1 = iterator->second.sizex1 / 2 - th;
536  G4double size_x2 = iterator->second.sizex2 / 2 - th;
537  G4double size_y1 = iterator->second.sizey1 / 2 - th;
538  G4double size_y2 = iterator->second.sizey2 / 2 - th;
539  G4double size_z = iterator->second.sizez / 2 - th;
540 
541  G4Trap* block_glass = new G4Trap("solid_glass",
542  size_z, // G4double pDz,
543  iterator->second.pTheta, 0, // G4double pTheta, G4double pPhi,
544  size_y1, size_x1, size_x1, // G4double pDy1, G4double pDx1, G4double pDx2,
545  0, // G4double pAlp1,
546  size_y2, size_x2, size_x2, // G4double pDy2, G4double pDx3, G4double pDx4,
547  0 // G4double pAlp2 //
548  );
549 
550  return block_glass;
551 }
552 
553 G4Trap* PHG4BarrelEcalDetector::GetTowerTrap(std::map<std::string, towerposition>::iterator iterator)
554 {
555  G4Trap* block_tower = new G4Trap("solid_tower",
556  iterator->second.sizez / 2, // G4double pDz,
557  iterator->second.pTheta, 0, // G4double pTheta, G4double pPhi,
558  iterator->second.sizey1 / 2, iterator->second.sizex1 / 2, iterator->second.sizex1 / 2., // G4double pDy1, G4double pDx1, G4double pDx2,
559  0, // G4double pAlp1,
560  iterator->second.sizey2 / 2, iterator->second.sizex2 / 2, iterator->second.sizex2 / 2., // G4double pDy2, G4double pDx3, G4double pDx4,
561  0 // G4double pAlp2 //
562  );
563 
564  return block_tower;
565 }
566 
567 G4Trap* PHG4BarrelEcalDetector::GetGlassTrapSubtract(std::map<std::string, towerposition>::iterator iterator)
568 {
569  G4double th = m_Params->get_double_param("thickness_wall") * cm;
570  G4double size_x1 = iterator->second.sizex1 / 2 - th + overlap;
571  G4double size_x2 = iterator->second.sizex2 / 2 - th + overlap;
572  G4double size_y1 = iterator->second.sizey1 / 2 - th + overlap;
573  G4double size_y2 = iterator->second.sizey2 / 2 - th + overlap;
574  G4double size_z = iterator->second.sizez / 2 - th / 2 + overlap;
575 
576  G4Trap* block_glass = new G4Trap("solid_glass",
577  size_z, // G4double pDz,
578  iterator->second.pTheta, 0, // G4double pTheta, G4double pPhi,
579  size_y1, size_x1, size_x1, // G4double pDy1, G4double pDx1, G4double pDx2,
580  0, // G4double pAlp1,
581  size_y2, size_x2, size_x2, // G4double pDy2, G4double pDx3, G4double pDx4,
582  0 // G4double pAlp2 //
583  );
584 
585  return block_glass;
586 }
587 
588 G4Trap* PHG4BarrelEcalDetector::GetSiTrap(std::map<std::string, towerposition>::iterator iterator)
589 {
590  G4double size_x1 = iterator->second.sizex2 / 2 - overlap;
591  G4double size_x2 = iterator->second.sizex2 / 2 - overlap;
592  G4double size_y1 = iterator->second.sizey2 / 2 - overlap;
593  G4double size_y2 = iterator->second.sizey2 / 2 - overlap;
594  G4double size_z = silicon_width_half;
595 
596  G4Trap* block_si = new G4Trap("Si",
597  size_z, // G4double pDz,
598  iterator->second.pTheta, 0, // G4double pTheta, G4double pPhi,
599  size_y1, size_x1, size_x1, // G4double pDy1, G4double pDx1, G4double pDx2,
600  0, // G4double pAlp1,
601  size_y2, size_x2, size_x2, // G4double pDy2, G4double pDx3, G4double pDx4,
602  0 // G4double pAlp2 //
603  );
604 
605  return block_si;
606 }
607 
608 G4LogicalVolume*
609 PHG4BarrelEcalDetector::ConstructSi(std::map<std::string, towerposition>::iterator iterator)
610 {
611  G4Trap* block_solid = GetSiTrap(iterator);
612  G4Material* material_si = GetDetectorMaterial("G4_POLYSTYRENE");
613  assert(material_si);
614  G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, material_si,
615  G4String(string(iterator->first) + string("_solid_Si")), 0, 0,
616  nullptr);
617  return block_logic;
618 }
619 
620 G4Trap* PHG4BarrelEcalDetector::GetKaptonTrap(std::map<std::string, towerposition>::iterator iterator)
621 {
622  G4double size_x1 = iterator->second.sizex2 / 2 - overlap;
623  G4double size_x2 = iterator->second.sizex2 / 2 - overlap;
624  G4double size_y1 = iterator->second.sizey2 / 2 - overlap;
625  G4double size_y2 = iterator->second.sizey2 / 2 - overlap;
626  G4double size_z = kapton_width_half;
627  G4Trap* block_si = new G4Trap("Kapton",
628  size_z, // G4double pDz,
629  iterator->second.pTheta, 0, // G4double pTheta, G4double pPhi,
630  size_y1, size_x1, size_x1, // G4double pDy1, G4double pDx1, G4double pDx2,
631  0, // G4double pAlp1,
632  size_y2, size_x2, size_x2, // G4double pDy2, G4double pDx3, G4double pDx4,
633  0 // G4double pAlp2 //
634  );
635  return block_si;
636 }
637 
638 G4LogicalVolume*
639 PHG4BarrelEcalDetector::ConstructKapton(std::map<std::string, towerposition>::iterator iterator)
640 {
641  G4Trap* block_solid = GetKaptonTrap(iterator);
642  G4Material* material_kapton = GetDetectorMaterial("G4_KAPTON");
643  assert(material_kapton);
644  G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, material_kapton,
645  G4String(string(iterator->first) + string("_solid_Si")), 0, 0,
646  nullptr);
647  return block_logic;
648 }
649 
650 G4Trap* PHG4BarrelEcalDetector::GetSIO2Trap(std::map<std::string, towerposition>::iterator iterator)
651 {
652  G4double size_x1 = iterator->second.sizex2 / 2 - overlap;
653  G4double size_x2 = iterator->second.sizex2 / 2 - overlap;
654  G4double size_y1 = iterator->second.sizey2 / 2 - overlap;
655  G4double size_y2 = iterator->second.sizey2 / 2 - overlap;
656  G4double size_z = SIO2_width_half;
657 
658  G4Trap* block_SIO2 = new G4Trap("SIO2",
659  size_z, // G4double pDz,
660  iterator->second.pTheta, 0, // G4double pTheta, G4double pPhi,
661  size_y1, size_x1, size_x1, // G4double pDy1, G4double pDx1, G4double pDx2,
662  0, // G4double pAlp1,
663  size_y2, size_x2, size_x2, // G4double pDy2, G4double pDx3, G4double pDx4,
664  0 // G4double pAlp2 //
665  );
666  return block_SIO2;
667 }
668 
669 G4LogicalVolume*
670 PHG4BarrelEcalDetector::ConstructSIO2(std::map<std::string, towerposition>::iterator iterator)
671 {
672  G4Trap* block_solid = GetSIO2Trap(iterator);
673  G4Material* material_SIO2 = GetDetectorMaterial("Quartz");
674  assert(material_SIO2);
675  G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, material_SIO2,
676  G4String(string(iterator->first) + string("_solid_material_SIO2")), 0, 0,
677  nullptr);
678  return block_logic;
679 }
680 
681 G4Trap* PHG4BarrelEcalDetector::GetCarbonTrap(std::map<std::string, towerposition>::iterator iterator)
682 {
683  G4double size_x1 = iterator->second.sizex2 / 2 - overlap;
684  G4double size_x2 = iterator->second.sizex2 / 2 - overlap;
685  G4double size_y1 = iterator->second.sizey2 / 2 - overlap;
686  G4double size_y2 = iterator->second.sizey2 / 2 - overlap;
687  G4double size_z = Carbon_width_half;
688 
689  G4Trap* block_SIO2 = new G4Trap("C",
690  size_z, // G4double pDz,
691  iterator->second.pTheta, 0, // G4double pTheta, G4double pPhi,
692  size_y1, size_x1, size_x1, // G4double pDy1, G4double pDx1, G4double pDx2,
693  0, // G4double pAlp1,
694  size_y2, size_x2, size_x2, // G4double pDy2, G4double pDx3, G4double pDx4,
695  0 // G4double pAlp2 //
696  );
697 
698  return block_SIO2;
699 }
700 
701 G4LogicalVolume*
702 PHG4BarrelEcalDetector::ConstructCarbon(std::map<std::string, towerposition>::iterator iterator)
703 {
704  G4Trap* block_solid = GetCarbonTrap(iterator);
705  G4Material* material_Carbon = GetDetectorMaterial("G4_C");
706  assert(material_Carbon);
707  G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, material_Carbon,
708  G4String(string(iterator->first) + string("_solid_material_C")), 0, 0,
709  nullptr);
710  return block_logic;
711 }