ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4HybridHomogeneousCalorimeterDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4HybridHomogeneousCalorimeterDetector.cc
1 #include "PHG4HybridHomogeneousCalorimeterDetector.h"
2 #include "PHG4HybridHomogeneousCalorimeterDisplayAction.h"
3 
4 #include <phparameter/PHParameters.h>
5 
6 #include <g4main/PHG4Detector.h> // for PHG4Detector
7 #include <g4main/PHG4DisplayAction.h> // for PHG4DisplayAction
8 #include <g4main/PHG4Subsystem.h>
9 
10 #include <phool/recoConsts.h>
11 
12 #include <Geant4/G4Box.hh>
13 #include <Geant4/G4Cons.hh>
14 #include <Geant4/G4Element.hh> // for G4Element
15 #include <Geant4/G4LogicalBorderSurface.hh>
16 #include <Geant4/G4LogicalSkinSurface.hh>
17 #include <Geant4/G4LogicalVolume.hh>
18 #include <Geant4/G4Material.hh>
19 #include <Geant4/G4MaterialPropertiesTable.hh>
20 #include <Geant4/G4OpticalSurface.hh>
21 #include <Geant4/G4PVPlacement.hh>
22 #include <Geant4/G4Polyhedra.hh> // for G4RotationMatrix
23 #include <Geant4/G4RotationMatrix.hh> // for G4RotationMatrix
24 #include <Geant4/G4String.hh> // for G4String
25 #include <Geant4/G4SubtractionSolid.hh>
26 #include <Geant4/G4SystemOfUnits.hh>
27 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
28 #include <Geant4/G4Transform3D.hh> // for G4Transform3D
29 #include <Geant4/G4Types.hh> // for G4double
30 #include <Geant4/G4UnionSolid.hh>
31 #include <Geant4/G4VisAttributes.hh>
32 
33 #include <TSystem.h>
34 
35 #include <cmath>
36 #include <cstdlib>
37 #include <iostream>
38 #include <sstream>
39 #include <utility> // for pair, make_pair
40 
41 class G4VSolid;
42 class PHCompositeNode;
43 
44 using namespace std;
45 
46 //_______________________________________________________________________
47 PHG4HybridHomogeneousCalorimeterDetector::PHG4HybridHomogeneousCalorimeterDetector(PHG4Subsystem* subsys, PHCompositeNode* Node, PHParameters* parameters, const std::string& dnam)
48  : PHG4Detector(subsys, Node, dnam)
49  , m_SuperDetector("NONE")
50  , m_Params(parameters)
51  , m_DisplayAction(dynamic_cast<PHG4HybridHomogeneousCalorimeterDisplayAction*>(subsys->GetDisplayAction()))
52  , _towerlogicnameprefix("HybridHomogeneousCalorimeterTower")
53  , m_IsActive(m_Params->get_int_param("active"))
54  , m_AbsorberActive(m_Params->get_int_param("absorberactive"))
55  , m_doLightProp(false)
56 {
57 }
58 
59 //_______________________________________________________________________
60 int PHG4HybridHomogeneousCalorimeterDetector::IsInCrystalCalorimeter(G4VPhysicalVolume* volume) const
61 {
62  if (m_IsActive)
63  {
64  if (m_ActiveVolumeSet.find(volume) != m_ActiveVolumeSet.end())
65  {
66  return GetCaloType();
67  }
68  }
69  if (m_AbsorberActive)
70  {
71  if (m_PassiveVolumeSet.find(volume) != m_PassiveVolumeSet.end())
72  {
73  return -1;
74  }
75  }
76  return 0;
77 }
78 
79 //_______________________________________________________________________
80 void PHG4HybridHomogeneousCalorimeterDetector::ConstructMe(G4LogicalVolume* logicWorld)
81 {
82  if (Verbosity() > 0)
83  {
84  cout << "PHG4HybridHomogeneousCalorimeterDetector: Begin Construction" << endl;
85  }
86 
87  if (m_Params->get_string_param("mappingtower").empty())
88  {
89  cout << "ERROR in PHG4HybridHomogeneousCalorimeterDetector: No tower mapping file specified. Abort detector construction." << endl;
90  cout << "Please run set_string_param(\"mappingtower\", std::string filename ) first." << endl;
91  exit(1);
92  }
93 
94  /* Read parameters for detector construction and mapping from file */
95  ParseParametersFromTable();
96 
97  // G4LogicalVolume* support_frame = ConstructSupportFrame(logicWorld);
98 
99  /* Construct single calorimeter tower */
100  G4LogicalVolume* singletower = ConstructTower();
101 
102  /* Place calorimeter tower within envelope */
103  // PlaceTower(eemc_envelope_log, singletower);
104  PlaceTower(logicWorld, singletower); //,support_frame);
105 
106  return;
107 }
108 //_______________________________________________________________________
109 G4LogicalVolume* PHG4HybridHomogeneousCalorimeterDetector::ConstructSupportFrame(G4LogicalVolume* eemcenvelope)
110 {
111  G4double frame_r_in = m_Params->get_double_param("rMin2") * cm;
112  G4double frame_r_out = m_Params->get_double_param("rMax2") * cm;
113  G4double frame_dz = m_Params->get_double_param("dz") * cm;
114  G4double frame_z_pos = m_Params->get_double_param("place_z") * cm;
115 
116  G4int nSides = 12;
117  const G4int nZSlices = 4;
118  G4double zPosSlice[nZSlices] = {0, 0.001, frame_dz - 0.001, frame_dz};
119  G4double rinPosSlice[nZSlices] = {frame_r_in, frame_r_in, frame_r_in, frame_r_in};
120  G4double routPosSlice[nZSlices] = {frame_r_out, frame_r_out, frame_r_out, frame_r_out};
121 
122  G4VSolid* frame_full_solid = new G4Polyhedra(G4String("frame_full_solid"),
123  0,
124  2 * M_PI * rad,
125  nSides,
126  nZSlices,
127  zPosSlice,
128  rinPosSlice,
129  routPosSlice);
130  G4LogicalVolume* frame_full_logic = new G4LogicalVolume(frame_full_solid, GetDetectorMaterial("G4_Fe"), "frame_full_logic", 0, 0, 0);
131  GetDisplayAction()->AddVolume(frame_full_logic, "WIP");
132 
133  G4RotationMatrix* rotFrame = new G4RotationMatrix();
134  rotFrame->rotateZ(M_PI / 12);
135  new G4PVPlacement(rotFrame, G4ThreeVector(0, 0, frame_z_pos - frame_dz / 2),
136  frame_full_logic,
137  "frame_full_placed",
138  eemcenvelope,
139  0, 0, OverlapCheck());
140 
141  return frame_full_logic;
142 }
143 
144 //_______________________________________________________________________
145 G4LogicalVolume* PHG4HybridHomogeneousCalorimeterDetector::ConstructTower()
146 {
147  if (Verbosity() > 0)
148  {
149  cout << "PHG4HybridHomogeneousCalorimeterDetector: Build logical volume for single tower..." << endl;
150  }
151 
152  /* dimensions of full tower */
153  G4double crystal_dx = m_Params->get_double_param("crystal_dx") * cm;
154  G4double crystal_dy = m_Params->get_double_param("crystal_dy") * cm;
155  G4double crystal_dz = m_Params->get_double_param("crystal_dz") * cm;
156 
157  G4double carbon_thickness = m_Params->get_double_param("carbon_gap") * cm;
158  G4double airgap_crystal_carbon = m_Params->get_double_param("air_gap") * cm;
159 
160  G4double reflective_foil_thickness = m_Params->get_double_param("reflective_foil_thickness") * cm;
161  G4double tedlar_thickness = m_Params->get_double_param("tedlar_thickness") * cm;
162  bool doWrapping = false;
163  if (reflective_foil_thickness > 0 || tedlar_thickness > 0) doWrapping = true;
164 
165  G4int sensor_count = m_Params->get_int_param("sensor_count");
166  G4double sensor_dimension = m_Params->get_double_param("sensor_dimension") * cm;
167  G4double sensor_thickness = m_Params->get_double_param("sensor_thickness") * cm;
168  bool doSensors = false;
169  if (sensor_dimension > 0 && sensor_count > 0 && sensor_thickness > 0) doSensors = true;
170 
171  G4double tower_dx = crystal_dx + 2 * (carbon_thickness + airgap_crystal_carbon + reflective_foil_thickness + tedlar_thickness);
172  G4double tower_dy = crystal_dy + 2 * (carbon_thickness + airgap_crystal_carbon + reflective_foil_thickness + tedlar_thickness);
173  G4double tower_dz = crystal_dz + 2 * (carbon_thickness);
174  if (doSensors) tower_dz = crystal_dz + 2 * (carbon_thickness) + sensor_thickness;
175  int carbon_frame_style = m_Params->get_int_param("carbon_frame_style");
176 
177  recoConsts* rc = recoConsts::instance();
178  G4Material* WorldMaterial = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
179 
180  /* create logical volume for single tower */
181  // Building the single tower mother volume first
182  // Then the crystal/sci-glass will be put in this volume
183  // The shell will also placed in. The rest space leave for the air gap
184 
185  G4VSolid* single_tower_solid = new G4Box(G4String("single_tower_solid"), tower_dx / 2.0, tower_dy / 2.0, tower_dz / 2.0);
186  G4LogicalVolume* single_tower_logic = new G4LogicalVolume(single_tower_solid, WorldMaterial, "single_tower_logic", 0, 0, 0);
187  GetDisplayAction()->AddVolume(single_tower_logic, "Invisible");
188  /* create geometry volume for crystal inside single_tower */
189  G4VSolid* solid_crystal = new G4Box(G4String("single_crystal_solid"),
190  crystal_dx / 2.0,
191  crystal_dy / 2.0,
192  crystal_dz / 2.0);
193 
194  G4Material* material_shell = GetCarbonFiber();
195  if (carbon_frame_style == 0 && carbon_thickness > 0)
196  {
197  /* create geometry volume for frame (carbon fiber shell) inside single_tower */
198  G4VSolid* Carbon_hunk_solid = new G4Box(G4String("Carbon_hunk_solid"),
199  tower_dx / 2.0,
200  tower_dy / 2.0,
201  ((tower_dz / 2.0) - 1 * mm));
202 
203  G4VSolid* cutout_solid = new G4Box(G4String("lead_solid"), crystal_dx / 2.0, crystal_dy / 2.0, crystal_dz);
204 
205  G4SubtractionSolid* Carbon_shell_solid = new G4SubtractionSolid(G4String("Carbon_Shell_solid"),
206  Carbon_hunk_solid,
207  cutout_solid,
208  0,
209  G4ThreeVector(0.00 * mm, 0.00 * mm, 0.00 * mm));
210 
211  G4LogicalVolume* Carbon_shell_logic = new G4LogicalVolume(Carbon_shell_solid, material_shell, "Carbon_shell_logic", 0, 0, 0);
212  GetDisplayAction()->AddVolume(Carbon_shell_logic, "CarbonShell");
213 
214  /* Place structural frame in logical tower volume */
215  // ostringstream name_shell;
216  // name_shell.str("");
217  // name_shell << _towerlogicnameprefix << "_single_absorber";
218  string name_shell = _towerlogicnameprefix + "_single_shell";
219  G4VPhysicalVolume* physvol_carbon = new G4PVPlacement(0, G4ThreeVector(0, 0, sensor_thickness / 2), Carbon_shell_logic, name_shell, single_tower_logic, 0, 0, OverlapCheck());
220  m_PassiveVolumeSet.insert(physvol_carbon);
221  }
222  else if (carbon_frame_style == 1)
223  {
224  // the carbon spacers should go a few cm deep between the crystals
225  G4double carbon_frame_depth = m_Params->get_double_param("carbon_frame_depth") * cm;
226  G4VSolid* Carbon_hunk_solid = new G4Box(G4String("Carbon_hunk_solid"),
227  tower_dx / 2.0,
228  tower_dy / 2.0,
229  (carbon_frame_depth / 2.0));
230  G4VSolid* cutout_solid = new G4Box(G4String("cutout_solid"), (tower_dx - 2 * (carbon_thickness)) / 2.0, (tower_dy - 2 * (carbon_thickness)) / 2.0, crystal_dz);
231 
232  G4SubtractionSolid* Carbon_gap_solid = new G4SubtractionSolid(G4String("Carbon_gap_solid"),
233  Carbon_hunk_solid,
234  cutout_solid,
235  0,
236  G4ThreeVector(0.00 * mm, 0.00 * mm, 0.00 * mm));
237  G4LogicalVolume* logic_gap = new G4LogicalVolume(Carbon_gap_solid, material_shell, "Carbon_gap_logic", 0, 0, 0);
238  GetDisplayAction()->AddVolume(logic_gap, "CarbonShell");
239 
240  string name_gap = _towerlogicnameprefix + "_single_shell";
241  G4VPhysicalVolume* physvol_carbon_gap_front = new G4PVPlacement(0, G4ThreeVector(0, 0, tower_dz / 2 - carbon_frame_depth / 2 - carbon_thickness), logic_gap, name_gap + "_gap_front", single_tower_logic, 0, 0, OverlapCheck());
242  m_PassiveVolumeSet.insert(physvol_carbon_gap_front);
243  G4VPhysicalVolume* physvol_carbon_gap_back = new G4PVPlacement(0, G4ThreeVector(0, 0, -tower_dz / 2 + carbon_frame_depth / 2 + carbon_thickness + sensor_thickness), logic_gap, name_gap + "_gap_back", single_tower_logic, 0, 0, OverlapCheck());
244  m_PassiveVolumeSet.insert(physvol_carbon_gap_back);
245 
246  G4double carbon_face_lip = m_Params->get_double_param("carbon_face_lip") * cm;
247 
248  G4VSolid* Carbon_hunk_solid_face = new G4Box(G4String("Carbon_hunk_solid_face"),
249  tower_dx / 2.0,
250  tower_dy / 2.0,
251  (carbon_thickness / 2.0));
252  G4VSolid* cutout_solid_face = new G4Box(G4String("cutout_solid_face"), (tower_dx - 2.0 * carbon_face_lip) / 2.0, (tower_dy - 2.0 * carbon_face_lip) / 2.0, (tower_dz));
253 
254  G4SubtractionSolid* Carbon_face_solid = new G4SubtractionSolid(G4String("Carbon_face_solid"),
255  Carbon_hunk_solid_face,
256  cutout_solid_face,
257  0,
258  G4ThreeVector(0.00 * mm, 0.00 * mm, 0.00 * mm));
259  G4LogicalVolume* logic_face = new G4LogicalVolume(Carbon_face_solid, material_shell, "Carbon_face_logic", 0, 0, 0);
260  GetDisplayAction()->AddVolume(logic_face, "CarbonShell");
261 
262  string name_face = _towerlogicnameprefix + "_single_shell";
263  G4VPhysicalVolume* physvol_carbon_face_front = new G4PVPlacement(0, G4ThreeVector(0, 0, tower_dz / 2 - carbon_thickness / 2), logic_face, name_face + "_face_front", single_tower_logic, 0, 0, OverlapCheck());
264  m_PassiveVolumeSet.insert(physvol_carbon_face_front);
265  G4VPhysicalVolume* physvol_carbon_face_back = new G4PVPlacement(0, G4ThreeVector(0, 0, -tower_dz / 2 + carbon_thickness / 2 + sensor_thickness), logic_face, name_face + "_face_back", single_tower_logic, 0, 0, OverlapCheck());
266  m_PassiveVolumeSet.insert(physvol_carbon_face_back);
267  }
268  // wrap towers in VM2000 and Tedlar foils if requested
269  if (doWrapping)
270  {
271  G4VSolid* VM2000_hunk_solid = new G4Box(G4String("VM2000_hunk_solid"),
272  (crystal_dx + 2 * reflective_foil_thickness) / 2.0,
273  (crystal_dy + 2 * reflective_foil_thickness) / 2.0,
274  ((crystal_dz / 2.0) - 1 * mm));
275 
276  G4VSolid* cutout_solid_VM2000 = new G4Box(G4String("cutout_solid_VM2000"), crystal_dx / 2.0, crystal_dy / 2.0, crystal_dz);
277 
278  G4SubtractionSolid* VM2000_foil_solid = new G4SubtractionSolid(G4String("VM2000_foil_solid"),
279  VM2000_hunk_solid,
280  cutout_solid_VM2000,
281  0,
282  G4ThreeVector(0.00 * mm, 0.00 * mm, 0.00 * mm));
283 
284  G4LogicalVolume* VM2000_foil_logic = new G4LogicalVolume(VM2000_foil_solid, GetVM2000Material(), "VM2000_foil_logic", 0, 0, 0);
285  GetDisplayAction()->AddVolume(VM2000_foil_logic, "VM2000");
286 
287  string name_VM2000foil = _towerlogicnameprefix + "_single_foil_VM2000";
288  G4VPhysicalVolume* physvol_VM2000 = new G4PVPlacement(0, G4ThreeVector(0, 0, sensor_thickness / 2), VM2000_foil_logic, name_VM2000foil, single_tower_logic, 0, 0, OverlapCheck());
289  m_PassiveVolumeSet.insert(physvol_VM2000);
290 
291  /* create geometry volume for frame (carbon fiber shell) inside single_tower */
292  G4VSolid* Tedlar_hunk_solid = new G4Box(G4String("Tedlar_hunk_solid"),
293  (crystal_dx + 2 * reflective_foil_thickness + 2 * tedlar_thickness) / 2.0,
294  (crystal_dy + 2 * reflective_foil_thickness + 2 * tedlar_thickness) / 2.0,
295  ((crystal_dz / 2.0) - 1 * mm));
296 
297  G4VSolid* cutout_solid_Tedlar = new G4Box(G4String("cutout_solid_Tedlar"), (crystal_dx + 2 * reflective_foil_thickness) / 2.0, (crystal_dy + 2 * reflective_foil_thickness) / 2.0, crystal_dz);
298 
299  G4SubtractionSolid* Tedlar_foil_solid = new G4SubtractionSolid(G4String("Tedlar_foil_solid"),
300  Tedlar_hunk_solid,
301  cutout_solid_Tedlar,
302  0,
303  G4ThreeVector(0.00 * mm, 0.00 * mm, 0.00 * mm));
304 
305  G4LogicalVolume* Tedlar_foil_logic = new G4LogicalVolume(Tedlar_foil_solid, GetTedlarMaterial(), "Tedlar_foil_logic", 0, 0, 0);
306  GetDisplayAction()->AddVolume(Tedlar_foil_logic, "Tedlar");
307 
308  string name_Tedlarfoil = _towerlogicnameprefix + "_single_foil_Tedlar";
309  G4VPhysicalVolume* physvol_Tedlar = new G4PVPlacement(0, G4ThreeVector(0, 0, sensor_thickness / 2), Tedlar_foil_logic, name_Tedlarfoil, single_tower_logic, 0, 0, OverlapCheck());
310  m_PassiveVolumeSet.insert(physvol_Tedlar);
311  }
312  /* create logical volumes for crystal inside single_tower */
313  G4double M_para = m_Params->get_double_param("material");
314  // set default material for hom calo
315  G4Material* material_Scin = GetDetectorMaterial("G4_PbWO4");
316  if (M_para > 0) material_Scin = GetScintillatorMaterial(M_para);
317 
318  if (m_doLightProp)
319  {
320  //scintillation and optical properties for the crystal
321  CrystalTable(material_Scin);
322  }
323 
324  G4LogicalVolume* logic_crystal = new G4LogicalVolume(solid_crystal, material_Scin, "single_crystal_logic", 0, 0, 0);
325 
326  if (m_doLightProp)
327  {
328  //crystal optical surface
329  SurfaceTable(logic_crystal);
330  }
331  GetDisplayAction()->AddVolume(logic_crystal, "Crystal");
332 
333  /* Place crystal in logical tower volume */
334  string name_crystal = _towerlogicnameprefix + "_single_crystal";
335  G4VPhysicalVolume* physvol_crys = new G4PVPlacement(0, G4ThreeVector(0, 0, sensor_thickness / 2), logic_crystal, name_crystal, single_tower_logic, 0, 0, OverlapCheck());
336  m_ActiveVolumeSet.insert(physvol_crys);
337 
338  if (doSensors)
339  {
340  G4VSolid* single_sensor_solid = new G4Box("single_sensor_solid", sensor_dimension / 2., sensor_dimension / 2., sensor_thickness / 2.);
341  G4Material* material_Sensor = GetDetectorMaterial("G4_Si");
342 
343  G4LogicalVolume* single_sensor_logic = new G4LogicalVolume(single_sensor_solid, material_Sensor, "single_sensor_logic", 0, 0, 0);
344  GetDisplayAction()->AddVolume(single_sensor_logic, "Sensor");
345 
346  string name_sensor = _towerlogicnameprefix + "_single_sensor";
347 
348  G4VPhysicalVolume* physvol_sensor_0 = new G4PVPlacement(0, G4ThreeVector(0, 0, -tower_dz / 2 + carbon_thickness + sensor_thickness / 2), single_sensor_logic, name_sensor, single_tower_logic, 0, 0, OverlapCheck());
349  m_ActiveVolumeSet.insert(physvol_sensor_0);
350  if (m_doLightProp)
351  {
352  MakeBoundary(physvol_crys, physvol_sensor_0);
353  }
354  }
355  if (Verbosity() > 0)
356  cout << "PHG4HybridHomogeneousCalorimeterDetector: Building logical volume for single tower done." << endl;
357 
358  return single_tower_logic;
359 }
360 
361 int PHG4HybridHomogeneousCalorimeterDetector::PlaceTower(G4LogicalVolume* eemcenvelope, G4LogicalVolume* singletower) //, G4LogicalVolume* framelogic)
362 {
363  /* Loop over all tower positions in vector and place tower */
364  typedef std::map<std::string, towerposition>::iterator it_type;
365  for (it_type iterator = _map_tower.begin(); iterator != _map_tower.end(); ++iterator)
366  {
367  if (Verbosity() > 0)
368  {
369  cout << "PHG4HybridHomogeneousCalorimeterDetector: Place tower " << iterator->first
370  << " idx_j = " << iterator->second.idx_j << ", idx_k = " << iterator->second.idx_k
371  << " at x = " << iterator->second.x << " , y = " << iterator->second.y << " , z = " << iterator->second.z << endl;
372  }
373  int copyno = (iterator->second.idx_j << 16) + iterator->second.idx_k;
374  new G4PVPlacement(0, G4ThreeVector(iterator->second.x, iterator->second.y, iterator->second.z),
375  singletower,
376  iterator->first,
377  eemcenvelope,
378  0, copyno, OverlapCheck());
379  }
380 
381  return 0;
382 }
383 
384 //_______________________________________________________________________
385 G4Material*
386 PHG4HybridHomogeneousCalorimeterDetector::GetScintillatorMaterial(float setting)
387 {
388  G4Element* ele_O = new G4Element("Oxygen", "O", 8., 16.00 * g / mole);
389  G4Element* ele_Si = new G4Element("Silicon", "Si", 14., 28.09 * g / mole);
390  G4Element* ele_B = new G4Element("Boron", "B", 5., 10.811 * g / mole);
391  G4Element* ele_Na = new G4Element("Sodium", "Na", 11., 22.99 * g / mole);
392  G4Element* ele_Mg = new G4Element("Magnesium", "Mg", 12., 24.30 * g / mole);
393  G4Element* ele_Pb = new G4Element("Lead", "Pb", 82., 207.2 * g / mole);
394  G4Element* ele_Ba = new G4Element("Barium", "Ba", 56., 137.3 * g / mole);
395  G4Element* ele_Gd = new G4Element("Gadolinium", "Gd", 64., 157.3 * g / mole);
396 
397  G4Material* material_Scin = GetDetectorMaterial("G4_PbWO4");
398 
399  if ((setting > 0.) && (setting < 1.))
400  {
401  material_Scin = GetDetectorMaterial("G4_PbWO4");
402  // g4MatData.push_back(0.0333333*mm/MeV);
403  if (Verbosity()) cout << "Set G4_PbWO4..." << endl;
404  }
405  else if ((setting > 1.) && (setting < 2.))
406  {
407  material_Scin = GetDetectorMaterial("G4_GLASS_LEAD");
408  if (Verbosity()) cout << "Set G4_GLASS_LEAD..." << endl;
409  }
410  else if ((setting > 2.) && (setting < 3.))
411  {
412  material_Scin = GetDetectorMaterial("G4_BARIUM_SULFATE");
413  if (Verbosity()) cout << "Set G4_BARIUM_SULFATE..." << endl;
414  }
415  else if ((setting > 3.) && (setting < 4.))
416  {
417  material_Scin = GetDetectorMaterial("G4_CESIUM_IODIDE");
418  if (Verbosity()) cout << "Set G4_CESIUM_IODIDE..." << endl;
419  }
420  else if ((setting > 4.) && (setting < 5.))
421  {
422  material_Scin = new G4Material("material_Scin", 4.5 * g / cm3, 5);
423  material_Scin->AddElement(ele_Si, 21.9 * perCent);
424  material_Scin->AddElement(ele_B, 8.8 * perCent);
425  material_Scin->AddElement(ele_Na, 10.4 * perCent);
426  material_Scin->AddElement(ele_Mg, 6.5 * perCent);
427  material_Scin->AddElement(ele_O, 52.4 * perCent);
428 
429  if (Verbosity()) cout << "Set Sciglass..." << endl;
430  }
431  else if ((setting > 5.) && (setting < 6.))
432  {
433  material_Scin = new G4Material("material_Scin", 9.0 * g / cm3, 5);
434  material_Scin->AddElement(ele_Si, 21.9 * perCent);
435  material_Scin->AddElement(ele_B, 8.8 * perCent);
436  material_Scin->AddElement(ele_Na, 10.4 * perCent);
437  material_Scin->AddElement(ele_Mg, 6.5 * perCent);
438  material_Scin->AddElement(ele_O, 52.4 * perCent);
439 
440  if (Verbosity()) cout << "Set heavier Sciglass..." << endl;
441  }
442  else if ((setting > 6.) && (setting < 7.))
443  {
444  material_Scin = new G4Material("material_Scin", 4.5 * g / cm3, 3);
445  material_Scin->AddElement(ele_Si, 21.9 * perCent);
446  material_Scin->AddElement(ele_O, 52.4 * perCent);
447  material_Scin->AddElement(ele_Pb, 25.7 * perCent);
448 
449  if (Verbosity()) cout << "Set Sciglass contained lead..." << endl;
450  }
451  else if ((setting > 7.) && (setting < 8.))
452  {
453  material_Scin = new G4Material("material_Scin", 4.22 * g / cm3, 4);
454  material_Scin->AddElement(ele_O, 0.261);
455  material_Scin->AddElement(ele_Ba, 0.3875);
456  material_Scin->AddElement(ele_Si, 0.1369);
457  material_Scin->AddElement(ele_Gd, 0.2146);
458 
459  if (Verbosity()) cout << "Set Sciglass from Nathaly" << endl;
460  }
461  else if ((setting > 8.) && (setting < 9.))
462  {
463  material_Scin = new G4Material("material_Scin", 3.8 * g / cm3, 3);
464  material_Scin->AddElement(ele_O, 0.293);
465  material_Scin->AddElement(ele_Ba, 0.502);
466  material_Scin->AddElement(ele_Si, 0.205);
467 
468  if (Verbosity()) cout << "Set Sciglass from g4e" << endl;
469  }
470 
471  return material_Scin;
472 }
473 
474 //_____________________________________________________________________________
475 void PHG4HybridHomogeneousCalorimeterDetector::CrystalTable(G4Material* mat)
476 {
477  //scintillation and optical properties
478 
479  //already done
480  if (mat->GetMaterialPropertiesTable()) return;
481 
482  //G4cout << "OpTable::CrystalTable, " << mat->GetMaterialPropertiesTable() << G4endl;
483 
484  const G4int ntab = 2;
485  G4double scin_en[] = {2.9 * eV, 3. * eV}; // 420 nm (the range is 414 - 428 nm)
486  G4double scin_fast[] = {1., 1.};
487 
488  G4MaterialPropertiesTable* tab = new G4MaterialPropertiesTable();
489 
490  tab->AddProperty("FASTCOMPONENT", scin_en, scin_fast, ntab);
491  tab->AddConstProperty("FASTTIMECONSTANT", 6 * ns);
492  tab->AddConstProperty("SCINTILLATIONYIELD", 200 / MeV); // 200/MEV nominal 10
493  tab->AddConstProperty("RESOLUTIONSCALE", 1.);
494 
495  G4double opt_en[] = {1.551 * eV, 3.545 * eV}; // 350 - 800 nm
496  G4double opt_r[] = {2.4, 2.4};
497  G4double opt_abs[] = {200 * cm, 200 * cm};
498 
499  tab->AddProperty("RINDEX", opt_en, opt_r, ntab);
500  tab->AddProperty("ABSLENGTH", opt_en, opt_abs, ntab);
501 
502  mat->SetMaterialPropertiesTable(tab);
503 
504 } //CrystalTable
505 
506 //_____________________________________________________________________________
507 void PHG4HybridHomogeneousCalorimeterDetector::SurfaceTable(G4LogicalVolume* vol)
508 {
509  //crystal optical surface
510 
511  G4OpticalSurface* surface = new G4OpticalSurface("CrystalSurface", unified, polished, dielectric_metal);
512  //G4LogicalSkinSurface *csurf =
513  new G4LogicalSkinSurface("CrystalSurfaceL", vol, surface);
514 
515  //surface material
516  const G4int ntab = 2;
517  G4double opt_en[] = {1.551 * eV, 3.545 * eV}; // 350 - 800 nm
518  G4double reflectivity[] = {0.8, 0.8};
519  G4double efficiency[] = {0.9, 0.9};
520  G4MaterialPropertiesTable* surfmat = new G4MaterialPropertiesTable();
521  surfmat->AddProperty("REFLECTIVITY", opt_en, reflectivity, ntab);
522  surfmat->AddProperty("EFFICIENCY", opt_en, efficiency, ntab);
523  surface->SetMaterialPropertiesTable(surfmat);
524  //csurf->DumpInfo();
525 
526 } //SurfaceTable
527 
528 //_____________________________________________________________________________
529 void PHG4HybridHomogeneousCalorimeterDetector::MakeBoundary(G4VPhysicalVolume* crystal, G4VPhysicalVolume* opdet)
530 {
531  //optical boundary between the crystal and optical photons detector
532 
533  G4OpticalSurface* surf = new G4OpticalSurface("OpDetS");
534  //surf->SetType(dielectric_dielectric); // photons go to the detector, must have rindex defined
535  surf->SetType(dielectric_metal); // photon is absorbed when reaching the detector, no material rindex required
536  //surf->SetFinish(ground);
537  surf->SetFinish(polished);
538  //surf->SetModel(unified);
539  surf->SetModel(glisur);
540 
541  new G4LogicalBorderSurface("OpDetB", crystal, opdet, surf);
542 
543  const G4int ntab = 2;
544  G4double opt_en[] = {1.551 * eV, 3.545 * eV}; // 350 - 800 nm
545  //G4double reflectivity[] = {0., 0.};
546  G4double reflectivity[] = {0.1, 0.1};
547  //G4double reflectivity[] = {0.9, 0.9};
548  //G4double reflectivity[] = {1., 1.};
549  G4double efficiency[] = {1., 1.};
550 
551  G4MaterialPropertiesTable* surfmat = new G4MaterialPropertiesTable();
552  surfmat->AddProperty("REFLECTIVITY", opt_en, reflectivity, ntab);
553  surfmat->AddProperty("EFFICIENCY", opt_en, efficiency, ntab);
554  surf->SetMaterialPropertiesTable(surfmat);
555 
556 } //MakeBoundary
557 
558 G4Material* PHG4HybridHomogeneousCalorimeterDetector::GetTedlarMaterial()
559 {
560  static string matname = "HybridHomogeneousTedlar";
561  G4Material* tedlar = GetDetectorMaterial(matname, false); // false suppresses warning that material does not exist
562  if (!tedlar)
563  {
564  G4double density_tedlar = 1.43 * g / cm3;
565  tedlar = new G4Material(matname, density_tedlar, 3);
566  tedlar->AddElement(G4Element::GetElement("C"), 2);
567  tedlar->AddElement(G4Element::GetElement("F"), 2);
568  tedlar->AddElement(G4Element::GetElement("H"), 2);
569  }
570  return tedlar;
571 }
572 
573 G4Material* PHG4HybridHomogeneousCalorimeterDetector::GetVM2000Material()
574 {
575  static string matname = "HybridHomogeneousVM2000";
576  G4Material* VM2000 = GetDetectorMaterial(matname, false); // false suppresses warning that material does not exist
577  if (!VM2000)
578  {
579  G4double density_VM2000 = 1.43 * g / cm3;
580  VM2000 = new G4Material(matname, density_VM2000, 3);
581  VM2000->AddElement(G4Element::GetElement("C"), 2);
582  VM2000->AddElement(G4Element::GetElement("F"), 2);
583  VM2000->AddElement(G4Element::GetElement("H"), 2);
584 
585  G4MaterialPropertiesTable* mptVM2000 = new G4MaterialPropertiesTable();
586  const G4int nEntriesVM2000 = 31;
587 
588  G4double photonE_VM2000[nEntriesVM2000] =
589  {1.37760 * eV, 1.45864 * eV, 1.54980 * eV, 1.65312 * eV, 1.71013 * eV, 1.77120 * eV, 1.83680 * eV, 1.90745 * eV, 1.98375 * eV, 2.06640 * eV, 2.10143 * eV, 2.13766 * eV, 2.17516 * eV, 2.21400 * eV, 2.25426 * eV, 2.29600 * eV, 2.33932 * eV, 2.38431 * eV, 2.43106 * eV, 2.47968 * eV, 2.53029 * eV, 2.58300 * eV, 2.63796 * eV, 2.69531 * eV, 2.75520 * eV, 2.81782 * eV, 2.88335 * eV, 2.95200 * eV, 3.09960 * eV, 3.54241 * eV, 4.13281 * eV};
590  G4double refractiveIndex_VM2000[nEntriesVM2000] =
591  {1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42};
592  mptVM2000->AddProperty("RINDEX", photonE_VM2000, refractiveIndex_VM2000, nEntriesVM2000);
593 
594  VM2000->SetMaterialPropertiesTable(mptVM2000);
595  }
596  return VM2000;
597 }
598 // <material name="PolyvinylChloride">
599 // <D value="1.3" unit="g/cm3"/>
600 // <fraction n="0.04838" ref="H"/>
601 // <fraction n="0.38436" ref="C"/>
602 // <fraction n="0.56726" ref="Cl"/>
603 
604 int PHG4HybridHomogeneousCalorimeterDetector::ParseParametersFromTable()
605 {
606  /* Open the datafile, if it won't open return an error */
607  ifstream istream_mapping(m_Params->get_string_param("mappingtower"));
608  if (!istream_mapping.is_open())
609  {
610  cout << "ERROR in PHG4HybridHomogeneousCalorimeterDetector: Failed to open mapping file " << m_Params->get_string_param("mappingtower") << endl;
611  gSystem->Exit(1);
612  }
613 
614  /* loop over lines in file */
615  string line_mapping;
616  while (getline(istream_mapping, line_mapping))
617  {
618  /* Skip lines starting with / including a '#' */
619  if (line_mapping.find("#") != string::npos)
620  {
621  if (Verbosity() > 0)
622  {
623  cout << "PHG4HybridHomogeneousCalorimeterDetector: SKIPPING line in mapping file: " << line_mapping << endl;
624  }
625  continue;
626  }
627 
628  istringstream iss(line_mapping);
629 
630  /* If line starts with keyword Tower, add to tower positions */
631  if (line_mapping.find("Tower ") != string::npos)
632  {
633  unsigned idx_j, idx_k, idx_l;
634  G4double pos_x, pos_y, pos_z;
635  G4double size_x, size_y, size_z;
636  G4double rot_x, rot_y, rot_z;
637  G4double dummy;
638  string dummys;
639 
640  /* read string- break if error */
641  if (!(iss >> dummys >> dummy >> idx_j >> idx_k >> idx_l >> pos_x >> pos_y >> pos_z >> size_x >> size_y >> size_z >> rot_x >> rot_y >> rot_z))
642  {
643  cout << "ERROR in PHG4HybridHomogeneousCalorimeterDetector: Failed to read line in mapping file " << m_Params->get_string_param("mappingtower") << endl;
644  gSystem->Exit(1);
645  }
646 
647  /* Construct unique name for tower */
648  /* Mapping file uses cm, this class uses mm for length */
649  ostringstream towername;
650  towername.str("");
651  towername << _towerlogicnameprefix << "_j_" << idx_j << "_k_" << idx_k;
652 
653  /* Add Geant4 units */
654  pos_x = pos_x * cm;
655  pos_y = pos_y * cm;
656  pos_z = pos_z * cm;
657 
658  /* insert tower into tower map */
659  towerposition tower_new;
660  tower_new.x = pos_x;
661  tower_new.y = pos_y;
662  tower_new.z = pos_z;
663  tower_new.idx_j = idx_j;
664  tower_new.idx_k = idx_k;
665  _map_tower.insert(make_pair(towername.str(), tower_new));
666  }
667  else
668  {
669  /* If this line is not a comment and not a tower, save parameter as string / value. */
670  string parname;
671  G4double parval;
672 
673  /* read string- break if error */
674  if (!(iss >> parname >> parval))
675  {
676  cout << "ERROR in PHG4HybridHomogeneousCalorimeterDetector: Failed to read line in mapping file " << m_Params->get_string_param("mappingtower") << endl;
677  gSystem->Exit(1);
678  }
679 
680  _map_global_parameter.insert(make_pair(parname, parval));
681  }
682  }
683 
684  /* Update member variables for global parameters based on parsed parameter file */
685  std::map<string, G4double>::iterator parit;
686 
687  parit = _map_global_parameter.find("Gcrystal_dx");
688  if (parit != _map_global_parameter.end())
689  {
690  m_Params->set_double_param("crystal_dx", parit->second); // in cm
691  }
692 
693  parit = _map_global_parameter.find("Gcrystal_dy");
694  if (parit != _map_global_parameter.end())
695  {
696  m_Params->set_double_param("crystal_dy", parit->second); // in cm
697  }
698 
699  parit = _map_global_parameter.find("Gcrystal_dz");
700  if (parit != _map_global_parameter.end())
701  {
702  m_Params->set_double_param("crystal_dz", parit->second); // in cm
703  }
704 
705  parit = _map_global_parameter.find("Gcarbon_gap");
706  if (parit != _map_global_parameter.end())
707  {
708  m_Params->set_double_param("carbon_gap", parit->second); // in cm
709  }
710 
711  parit = _map_global_parameter.find("Gair_gap");
712  if (parit != _map_global_parameter.end())
713  {
714  m_Params->set_double_param("air_gap", parit->second); // in cm
715  }
716 
717  parit = _map_global_parameter.find("Gr1_inner");
718  if (parit != _map_global_parameter.end())
719  {
720  m_Params->set_double_param("rMin1", parit->second);
721  }
722 
723  parit = _map_global_parameter.find("Gr1_outer");
724  if (parit != _map_global_parameter.end())
725  {
726  m_Params->set_double_param("rMax1", parit->second);
727  }
728 
729  parit = _map_global_parameter.find("Gr2_inner");
730  if (parit != _map_global_parameter.end())
731  {
732  m_Params->set_double_param("rMin2", parit->second);
733  }
734 
735  parit = _map_global_parameter.find("Gr2_outer");
736  if (parit != _map_global_parameter.end())
737  {
738  m_Params->set_double_param("rMax2", parit->second);
739  }
740 
741  parit = _map_global_parameter.find("Gdz");
742  if (parit != _map_global_parameter.end())
743  {
744  m_Params->set_double_param("dz", parit->second);
745  }
746 
747  parit = _map_global_parameter.find("Gx0");
748  if (parit != _map_global_parameter.end())
749  {
750  m_Params->set_double_param("place_x", parit->second);
751  }
752 
753  parit = _map_global_parameter.find("Gy0");
754  if (parit != _map_global_parameter.end())
755  {
756  m_Params->set_double_param("place_y", parit->second);
757  }
758 
759  parit = _map_global_parameter.find("Gz0");
760  if (parit != _map_global_parameter.end())
761  {
762  m_Params->set_double_param("place_z", parit->second);
763  }
764 
765  parit = _map_global_parameter.find("Grot_x");
766  if (parit != _map_global_parameter.end())
767  {
768  m_Params->set_double_param("rot_x", parit->second * rad / deg);
769  }
770 
771  parit = _map_global_parameter.find("Grot_y");
772  if (parit != _map_global_parameter.end())
773  {
774  m_Params->set_double_param("rot_y", parit->second * rad / deg);
775  }
776 
777  parit = _map_global_parameter.find("Grot_z");
778  if (parit != _map_global_parameter.end())
779  {
780  m_Params->set_double_param("rot_z", parit->second * rad / deg);
781  }
782 
783  parit = _map_global_parameter.find("Gmaterial");
784  if (parit != _map_global_parameter.end())
785  {
786  m_Params->set_double_param("material", parit->second);
787  }
788 
789  parit = _map_global_parameter.find("Gcolor_R");
790  if (parit != _map_global_parameter.end())
791  {
792  m_Params->set_double_param("color_R", parit->second);
793  }
794 
795  parit = _map_global_parameter.find("Gcolor_G");
796  if (parit != _map_global_parameter.end())
797  {
798  m_Params->set_double_param("color_G", parit->second);
799  }
800 
801  parit = _map_global_parameter.find("Gcolor_B");
802  if (parit != _map_global_parameter.end())
803  {
804  m_Params->set_double_param("color_B", parit->second);
805  }
806  parit = _map_global_parameter.find("carbon_frame_style");
807  if (parit != _map_global_parameter.end())
808  {
809  m_Params->set_int_param("carbon_frame_style", parit->second);
810  }
811  parit = _map_global_parameter.find("carbon_frame_depth");
812  if (parit != _map_global_parameter.end())
813  {
814  m_Params->set_double_param("carbon_frame_depth", parit->second);
815  }
816  parit = _map_global_parameter.find("carbon_face_lip");
817  if (parit != _map_global_parameter.end())
818  {
819  m_Params->set_double_param("carbon_face_lip", parit->second);
820  }
821  parit = _map_global_parameter.find("Gtedlar_thickness");
822  if (parit != _map_global_parameter.end())
823  {
824  m_Params->set_double_param("tedlar_thickness", parit->second);
825  }
826  parit = _map_global_parameter.find("Greflective_foil_thickness");
827  if (parit != _map_global_parameter.end())
828  {
829  m_Params->set_double_param("reflective_foil_thickness", parit->second);
830  }
831  parit = _map_global_parameter.find("sensor_dimension");
832  if (parit != _map_global_parameter.end())
833  {
834  m_Params->set_double_param("sensor_dimension", parit->second);
835  }
836  parit = _map_global_parameter.find("sensor_count");
837  if (parit != _map_global_parameter.end())
838  {
839  m_Params->set_int_param("sensor_count", parit->second);
840  }
841  parit = _map_global_parameter.find("sensor_thickness");
842  if (parit != _map_global_parameter.end())
843  {
844  m_Params->set_double_param("sensor_thickness", parit->second);
845  }
846 
847  return 0;
848 }
849 
850 G4Material* PHG4HybridHomogeneousCalorimeterDetector::GetCarbonFiber()
851 {
852  static string matname = "HybridHomogeneousCarbonFiber";
853  G4Material* carbonfiber = GetDetectorMaterial(matname, false); // false suppresses warning that material does not exist
854  if (!carbonfiber)
855  {
856  G4double density_carbon_fiber = 1.44 * g / cm3;
857  carbonfiber = new G4Material(matname, density_carbon_fiber, 1);
858  carbonfiber->AddElement(G4Element::GetElement("C"), 1);
859  }
860  return carbonfiber;
861 }