ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4ForwardEcalDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4ForwardEcalDetector.cc
2 
4 
5 #include <phparameter/PHParameters.h>
6 
9 
10 #include <g4main/PHG4Detector.h> // for PHG4Detector
11 #include <g4main/PHG4DisplayAction.h> // for PHG4DisplayAction
12 #include <g4main/PHG4Subsystem.h>
13 
14 #include <phool/recoConsts.h>
15 
16 #include <Geant4/G4Box.hh>
17 #include <Geant4/G4Cons.hh>
18 #include <Geant4/G4LogicalSkinSurface.hh>
19 #include <Geant4/G4LogicalVolume.hh>
20 #include <Geant4/G4Material.hh>
21 #include <Geant4/G4MaterialPropertiesTable.hh>
22 #include <Geant4/G4OpticalSurface.hh>
23 #include <Geant4/G4PVPlacement.hh>
24 #include <Geant4/G4PVReplica.hh>
25 #include <Geant4/G4PhysicalConstants.hh>
26 #include <Geant4/G4Polyhedra.hh>
27 #include <Geant4/G4RotationMatrix.hh> // for G4RotationMatrix
28 #include <Geant4/G4String.hh> // for G4String
29 #include <Geant4/G4SubtractionSolid.hh>
30 #include <Geant4/G4SystemOfUnits.hh>
31 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
32 #include <Geant4/G4Transform3D.hh> // for G4Transform3D
33 #include <Geant4/G4Tubs.hh>
34 #include <Geant4/G4Types.hh> // for G4double, G4int
35 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
36 
37 #include <TSystem.h>
38 
39 #include <cassert>
40 #include <cmath>
41 #include <fstream>
42 #include <iostream>
43 #include <sstream>
44 #include <utility> // for pair, make_pair
45 
46 class G4VSolid;
47 class PHCompositeNode;
48 
49 //_______________________________________________________________________
51  : PHG4Detector(subsys, Node, dnam)
52  , m_DisplayAction(dynamic_cast<PHG4ForwardEcalDisplayAction*>(subsys->GetDisplayAction()))
53  , m_Params(parameters)
54  , m_GdmlConfig(PHG4GDMLUtility::GetOrMakeConfigNode(Node))
55  , m_ActiveFlag(m_Params->get_int_param("active"))
56  , m_AbsorberActiveFlag(m_Params->get_int_param("absorberactive"))
57  , m_doLightProp(false)
58 {
59  for (int i = 0; i < 3; i++)
60  {
61  m_TowerDx[i] = 30 * mm;
62  m_TowerDy[i] = 30 * mm;
63  m_TowerDz[i] = 170.0 * mm;
64  }
65  for (int i = 3; i < 7; i++)
66  {
67  m_TowerDx[i] = NAN;
68  m_TowerDy[i] = NAN;
69  m_TowerDz[i] = NAN;
70  }
71  m_RMin[0] = 110 * mm;
72  m_RMax[0] = 2250 * mm;
73  m_RMin[1] = 120 * mm;
74  m_RMax[1] = 2460 * mm;
75  m_Params->set_double_param("xoffset", 0.);
76  m_Params->set_double_param("yoffset", 0.);
77 
78  assert(m_GdmlConfig);
79 }
80 
81 //_______________________________________________________________________
83 {
84  G4LogicalVolume* mylogvol = volume->GetLogicalVolume();
85  if (m_ActiveFlag)
86  {
87  if (m_ScintiLogicalVolSet.find(mylogvol) != m_ScintiLogicalVolSet.end())
88  {
89  return 1;
90  }
91  }
93  {
94  if (m_AbsorberLogicalVolSet.find(mylogvol) != m_AbsorberLogicalVolSet.end())
95  {
96  return -1;
97  }
98  }
99  return 0;
100 }
101 
102 //_______________________________________________________________________
104 {
105  if (Verbosity() > 0)
106  {
107  std::cout << "PHG4ForwardEcalDetector: Begin Construction" << std::endl;
108  }
109 
110  /* Read parameters for detector construction and mappign from file */
112 
113  /* Create the cone envelope = 'world volume' for the crystal calorimeter */
116 
117  G4double tower_readout_dz = m_Params->get_double_param("tower_readout_dz") * cm;
118 
119  G4VSolid* beampipe_cutout = new G4Cons("FEMC_beampipe_cutout",
120  0, m_RMin[0],
121  0, m_RMin[1],
122  (m_dZ + tower_readout_dz),
123  0, 2 * M_PI);
124  G4VSolid* ecal_envelope_solid = new G4Cons("FEMC_envelope_solid_cutout",
125  0, m_RMax[0],
126  0, m_RMax[1],
127  (m_dZ + tower_readout_dz) / 2.0,
128  0, 2 * M_PI);
129  ecal_envelope_solid = new G4SubtractionSolid(G4String("hFEMC_envelope_solid"), ecal_envelope_solid, beampipe_cutout, 0, G4ThreeVector(m_Params->get_double_param("xoffset") * cm, m_Params->get_double_param("yoffset") * cm, 0.));
130 
131  G4LogicalVolume* ecal_envelope_log = new G4LogicalVolume(ecal_envelope_solid, WorldMaterial, "hFEMC_envelope", 0, 0, 0);
132 
133  /* Define visualization attributes for envelope cone */
134  GetDisplayAction()->AddVolume(ecal_envelope_log, "Envelope");
135 
136  /* Define rotation attributes for envelope cone */
137  G4RotationMatrix ecal_rotm;
138  ecal_rotm.rotateX(m_XRot);
139  ecal_rotm.rotateY(m_YRot);
140  ecal_rotm.rotateZ(m_ZRot);
141 
142  /* Place envelope cone in simulation */
143  std::string name_envelope = m_TowerLogicNamePrefix + "_envelope";
144 
145  new G4PVPlacement(G4Transform3D(ecal_rotm, G4ThreeVector(m_PlaceX, m_PlaceY, m_PlaceZ - tower_readout_dz / 2)),
146  ecal_envelope_log, name_envelope, logicWorld, 0, false, OverlapCheck());
147 
148  /* Construct single calorimeter towers */
149  G4LogicalVolume* singletower[7] = {nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr};
150  typedef std::map<std::string, towerposition>::iterator it_type;
151  for (it_type iterator = m_TowerPositionMap.begin(); iterator != m_TowerPositionMap.end(); ++iterator)
152  {
153  for (int i = 0; i < 7; i++)
154  {
155  if (iterator->second.type == i && singletower[i] == nullptr)
156  {
157  singletower[i] = ConstructTower(i);
158  }
159  }
160  }
161 
162  if (Verbosity() > 1)
163  {
164  std::cout << singletower << std::endl;
165  }
166  /* Place calorimeter towers within envelope */
167  PlaceTower(ecal_envelope_log, singletower);
168 
169  return;
170 }
171 
172 //_______________________________________________________________________
175 {
176  if (Verbosity() > 0)
177  {
178  std::cout << "PHG4ForwardEcalDetector: Build logical volume for single tower, type = " << type << std::endl;
179  }
180  /* create logical volume for single tower */
183 
184  /* create geometry volumes for scintillator and absorber plates to place inside single_tower */
185  // PHENIX EMCal JGL 3/27/2016
186  G4int nlayers = 66;
187  G4double thickness_layer = m_TowerDz[2] / (float) nlayers;
188  // update layer thickness with https://doi.org/10.1016/S0168-9002(02)01954-X
189  G4double thickness_cell = 5.6 * mm;
190  G4double thickness_absorber = 1.55 * mm; // 1.55mm absorber
191  G4double thickness_scintillator = 4.0 * mm; // 4mm scintillator
192  // notched in TiO2 coating in scintillator plate
193  G4double width_coating = m_Params->get_double_param("width_coating") * cm; // 4mm scintillator
194  G4Material* material_scintillator = GetScintillatorMaterial(); //GetDetectorMaterial("G4_POLYSTYRENE");
195  G4Material* material_absorber = GetDetectorMaterial("G4_Pb");
196  G4int nFibers = m_Params->get_int_param("nFibers");
197  G4double fiber_diam = m_Params->get_double_param("fiber_diam") * cm; // 4mm scintillator
198  G4double fiber_extra_length = 0.0;
199  // additional fiber length at end of calorimeter
200  if (fiber_diam > 0) fiber_extra_length = 0.5 * cm;
201  // width of plate in front of FEMC for clamping all layers
202  G4double clamp_plate_width = m_Params->get_double_param("clamp_plate_width") * cm;
203  // depth of readout (4cm)
204  G4double tower_readout_dz = m_Params->get_double_param("tower_readout_dz") * cm;
205 
206  G4VSolid* single_tower_solid = new G4Box("single_tower_solid2",
207  m_TowerDx[2] / 2.0,
208  m_TowerDy[2] / 2.0,
209  (m_TowerDz[2] + tower_readout_dz) / 2.0);
210  G4VSolid* single_tower_solid_replica = new G4Box("single_tower_solid_replica",
211  m_TowerDx[2] / 2.0,
212  m_TowerDy[2] / 2.0,
213  m_TowerDz[2] / 2.0);
214 
215  G4LogicalVolume* single_tower_logic = new G4LogicalVolume(single_tower_solid,
216  WorldMaterial,
217  "single_tower_logic2",
218  0, 0, 0);
219 
220  GetDisplayAction()->AddVolume(single_tower_logic, "SingleTower");
221 
222  if (Verbosity())
223  {
224  std::cout << " m_TowerDz[2] = " << m_TowerDz[2] << " thickness_layer = " << thickness_layer << " thickness_cell = " << thickness_cell << std::endl;
225  }
226 
227  if (thickness_layer <= thickness_cell)
228  {
229  std::cout << __PRETTY_FUNCTION__
230  << "Tower size z (m_TowerDz[2) from database is too thin. "
231  << "It does not fit the layer structure as described in https://doi.org/10.1016/S0168-9002(02)01954-X !" << std::endl
232  << "Abort" << std::endl;
233  std::cout << " m_TowerDz[2] = " << m_TowerDz[2] << " i.e. nlayers " << nlayers << " * thickness_layer " << thickness_layer << " <= thickness_cell " << thickness_cell << std::endl;
234  exit(1);
235  }
236 
237  //**********************************************************************************************
238  /* create logical and geometry volumes for minitower read-out unit */
239  //**********************************************************************************************
240  G4VSolid* miniblock_solid = new G4Box("miniblock_solid",
241  m_TowerDx[2] / 2.0,
242  m_TowerDy[2] / 2.0,
243  thickness_cell / 2.0);
244  G4LogicalVolume* miniblock_logic = new G4LogicalVolume(miniblock_solid,
245  WorldMaterial,
246  "miniblock_logic",
247  0, 0, 0);
248  GetDisplayAction()->AddVolume(miniblock_logic, "miniblock");
249  //**********************************************************************************************
250  /* create logical & geometry volumes for scintillator and absorber plates to place inside mini read-out unit */
251  //**********************************************************************************************
252  G4VSolid* solid_absorber = new G4Box("single_plate_absorber_solid2",
253  m_TowerDx[2] / 2.0,
254  m_TowerDy[2] / 2.0,
255  thickness_absorber / 2.0);
256 
257  G4VSolid* solid_scintillator = new G4Box("single_plate_scintillator2",
258  (m_TowerDx[2]) / 2.0,
259  (m_TowerDy[2]) / 2.0,
260  // (m_TowerDx[2] - 2*width_coating) / 2.0,
261  // (m_TowerDy[2] - 2*width_coating) / 2.0,
262  thickness_scintillator / 2.0);
263 
264  if (clamp_plate_width > 0)
265  {
266  G4VSolid* solid_clamp1 = new G4Box("single_plate_clamp_solid1",
267  m_TowerDx[2] / 2.0,
268  m_TowerDy[2] / 2.0,
269  clamp_plate_width / 2.0);
270  G4LogicalVolume* logic_clampplate = new G4LogicalVolume(solid_clamp1,
271  GetDetectorMaterial("G4_Fe"),
272  "logic_clampplate",
273  0, 0, 0);
274  m_AbsorberLogicalVolSet.insert(logic_clampplate);
275  GetDisplayAction()->AddVolume(logic_clampplate, "Clamp");
276  std::string name_clamp = m_TowerLogicNamePrefix + "_single_plate_clamp";
277 
278  new G4PVPlacement(0, G4ThreeVector(0, 0, (m_dZ + tower_readout_dz) / 2.0 - clamp_plate_width / 2.0),
279  logic_clampplate,
280  name_clamp,
281  single_tower_logic,
282  0, 0, OverlapCheck());
283  }
284  if (nFibers > 0 && nFibers == 5 && fiber_diam > 0)
285  {
286  G4VSolid* cutoutfiber_solid = new G4Tubs("cutoutfiber_solid",
287  0.0, 1.01 * fiber_diam / 2.0, m_TowerDz[2], 0.0, 2 * M_PI);
288 
289  single_tower_solid_replica = new G4SubtractionSolid(G4String("single_tower_solid_replica_cu1"), single_tower_solid_replica, cutoutfiber_solid, 0, G4ThreeVector(0, 0, 0.));
290  single_tower_solid_replica = new G4SubtractionSolid(G4String("single_tower_solid_replica_cu2"), single_tower_solid_replica, cutoutfiber_solid, 0, G4ThreeVector(-m_TowerDx[2] / 4.0, m_TowerDy[2] / 4.0, 0.));
291  single_tower_solid_replica = new G4SubtractionSolid(G4String("single_tower_solid_replica_cu3"), single_tower_solid_replica, cutoutfiber_solid, 0, G4ThreeVector(m_TowerDx[2] / 4.0, m_TowerDy[2] / 4.0, 0.));
292  single_tower_solid_replica = new G4SubtractionSolid(G4String("single_tower_solid_replica_cu4"), single_tower_solid_replica, cutoutfiber_solid, 0, G4ThreeVector(m_TowerDx[2] / 4.0, -m_TowerDy[2] / 4.0, 0.));
293  single_tower_solid_replica = new G4SubtractionSolid(G4String("single_tower_solid_replica_cu5"), single_tower_solid_replica, cutoutfiber_solid, 0, G4ThreeVector(-m_TowerDx[2] / 4.0, -m_TowerDy[2] / 4.0, 0.));
294 
295  solid_absorber = new G4SubtractionSolid(G4String("solid_absorber_cu1"), solid_absorber, cutoutfiber_solid, 0, G4ThreeVector(0, 0, 0.));
296  solid_absorber = new G4SubtractionSolid(G4String("solid_absorber_cu2"), solid_absorber, cutoutfiber_solid, 0, G4ThreeVector(-m_TowerDx[2] / 4.0, m_TowerDy[2] / 4.0, 0.));
297  solid_absorber = new G4SubtractionSolid(G4String("solid_absorber_cu3"), solid_absorber, cutoutfiber_solid, 0, G4ThreeVector(m_TowerDx[2] / 4.0, m_TowerDy[2] / 4.0, 0.));
298  solid_absorber = new G4SubtractionSolid(G4String("solid_absorber_cu4"), solid_absorber, cutoutfiber_solid, 0, G4ThreeVector(m_TowerDx[2] / 4.0, -m_TowerDy[2] / 4.0, 0.));
299  solid_absorber = new G4SubtractionSolid(G4String("solid_absorber_cu5"), solid_absorber, cutoutfiber_solid, 0, G4ThreeVector(-m_TowerDx[2] / 4.0, -m_TowerDy[2] / 4.0, 0.));
300 
301  solid_scintillator = new G4SubtractionSolid(G4String("solid_scintillator_cu1"), solid_scintillator, cutoutfiber_solid, 0, G4ThreeVector(0, 0, 0.));
302  solid_scintillator = new G4SubtractionSolid(G4String("solid_scintillator_cu2"), solid_scintillator, cutoutfiber_solid, 0, G4ThreeVector(-m_TowerDx[2] / 4.0, m_TowerDy[2] / 4.0, 0.));
303  solid_scintillator = new G4SubtractionSolid(G4String("solid_scintillator_cu3"), solid_scintillator, cutoutfiber_solid, 0, G4ThreeVector(m_TowerDx[2] / 4.0, m_TowerDy[2] / 4.0, 0.));
304  solid_scintillator = new G4SubtractionSolid(G4String("solid_scintillator_cu4"), solid_scintillator, cutoutfiber_solid, 0, G4ThreeVector(m_TowerDx[2] / 4.0, -m_TowerDy[2] / 4.0, 0.));
305  solid_scintillator = new G4SubtractionSolid(G4String("solid_scintillator_cu5"), solid_scintillator, cutoutfiber_solid, 0, G4ThreeVector(-m_TowerDx[2] / 4.0, -m_TowerDy[2] / 4.0, 0.));
306 
307  if (clamp_plate_width > 0)
308  {
309  G4VSolid* solid_clamp2 = new G4Box("single_plate_clamp_solid2",
310  m_TowerDx[2] / 2.0,
311  m_TowerDy[2] / 2.0,
312  clamp_plate_width / 2.0);
313  solid_clamp2 = new G4SubtractionSolid(G4String("solid_clamp2_cu1"), solid_clamp2, cutoutfiber_solid, 0, G4ThreeVector(0, 0, 0.));
314  solid_clamp2 = new G4SubtractionSolid(G4String("solid_clamp2_cu2"), solid_clamp2, cutoutfiber_solid, 0, G4ThreeVector(-m_TowerDx[2] / 4.0, m_TowerDy[2] / 4.0, 0.));
315  solid_clamp2 = new G4SubtractionSolid(G4String("solid_clamp2_cu3"), solid_clamp2, cutoutfiber_solid, 0, G4ThreeVector(m_TowerDx[2] / 4.0, m_TowerDy[2] / 4.0, 0.));
316  solid_clamp2 = new G4SubtractionSolid(G4String("solid_clamp2_cu4"), solid_clamp2, cutoutfiber_solid, 0, G4ThreeVector(m_TowerDx[2] / 4.0, -m_TowerDy[2] / 4.0, 0.));
317  solid_clamp2 = new G4SubtractionSolid(G4String("solid_clamp2_cu5"), solid_clamp2, cutoutfiber_solid, 0, G4ThreeVector(-m_TowerDx[2] / 4.0, -m_TowerDy[2] / 4.0, 0.));
318  G4LogicalVolume* logic_clampplate2 = new G4LogicalVolume(solid_clamp2,
319  GetDetectorMaterial("G4_Fe"),
320  "logic_clampplate2",
321  0, 0, 0);
322  m_AbsorberLogicalVolSet.insert(logic_clampplate2);
323  GetDisplayAction()->AddVolume(logic_clampplate2, "Clamp");
324  std::string name_clamp = m_TowerLogicNamePrefix + "_single_plate_clamp2";
325 
326  new G4PVPlacement(0, G4ThreeVector(0, 0, (tower_readout_dz) / 2.0 - clamp_plate_width - m_TowerDz[2] / 2.0 - clamp_plate_width / 2.0),
327  logic_clampplate2,
328  name_clamp,
329  single_tower_logic,
330  0, 0, OverlapCheck());
331  }
332  }
333 
334  if (width_coating > 0)
335  {
336  G4double depthCoating = 0.93 * thickness_scintillator;
337  G4double zPlaneCoating[2] = {0, depthCoating};
338  G4double rInnerCoating[2] = {(m_TowerDx[2] - 2 * width_coating) / 2, (m_TowerDx[2] - 2 * width_coating) / 2.0};
339  G4double rOuterCoating[2] = {m_TowerDx[2] / 2.0, m_TowerDx[2] / 2.0};
340  G4VSolid* solid_coating = new G4Polyhedra("solid_coating",
341  0, 2 * M_PI,
342  4, 2,
343  zPlaneCoating, rInnerCoating, rOuterCoating);
344  G4LogicalVolume* logic_coating = new G4LogicalVolume(solid_coating,
346  "logic_coating",
347  0, 0, 0);
348  m_AbsorberLogicalVolSet.insert(logic_coating);
349  GetDisplayAction()->AddVolume(logic_coating, "Coating");
350  std::string name_coating = m_TowerLogicNamePrefix + "_single_plate_coating";
351 
352  G4RotationMatrix* rotCoating = new G4RotationMatrix();
353  rotCoating->rotateZ(M_PI / 4);
354  new G4PVPlacement(rotCoating, G4ThreeVector(0, 0, thickness_cell / 2.0 - depthCoating),
355  logic_coating,
356  name_coating,
357  miniblock_logic,
358  0, 0, OverlapCheck());
359  solid_scintillator = new G4SubtractionSolid(G4String("solid_scintillator_cuCoating"), solid_scintillator, solid_coating, rotCoating, G4ThreeVector(0, 0, thickness_scintillator / 2.0 - depthCoating));
360  }
361  G4LogicalVolume* logic_absorber = new G4LogicalVolume(solid_absorber,
362  material_absorber,
363  "single_plate_absorber_logic2",
364  0, 0, 0);
365  m_AbsorberLogicalVolSet.insert(logic_absorber);
366  G4LogicalVolume* logic_scint = new G4LogicalVolume(solid_scintillator,
367  material_scintillator,
368  "hEcal_scintillator_plate_logic2",
369  0, 0, 0);
370  m_ScintiLogicalVolSet.insert(logic_scint);
371  if (m_doLightProp)
372  {
373  SurfaceTable(logic_scint);
374  }
375 
376  GetDisplayAction()->AddVolume(logic_absorber, "Absorber");
377  GetDisplayAction()->AddVolume(logic_scint, "Scintillator");
378 
379  std::string name_absorber = m_TowerLogicNamePrefix + "_single_plate_absorber2";
380  std::string name_scintillator = m_TowerLogicNamePrefix + "_single_plate_scintillator2";
381 
382  new G4PVPlacement(0, G4ThreeVector(0, 0, -thickness_cell / 2.0 + thickness_absorber / 2.0),
383  logic_absorber,
384  name_absorber,
385  miniblock_logic,
386  0, 0, OverlapCheck());
387 
388  new G4PVPlacement(0, G4ThreeVector(0, 0, thickness_cell / 2.0 - thickness_scintillator / 2.0),
389  logic_scint,
390  name_scintillator,
391  miniblock_logic,
392  0, 0, OverlapCheck());
393 
394  G4LogicalVolume* single_tower_replica_logic = new G4LogicalVolume(single_tower_solid_replica,
395  WorldMaterial,
396  "single_tower_replica_logic",
397  0, 0, 0);
398  /* create replica within tower */
399  std::string name_tower = m_TowerLogicNamePrefix;
400  new G4PVReplica(name_tower, miniblock_logic, single_tower_replica_logic,
401  kZAxis, nlayers, thickness_layer, 0);
402 
403  GetDisplayAction()->AddVolume(single_tower_replica_logic, "SingleTower");
404 
405  new G4PVPlacement(0, G4ThreeVector(0, 0, (tower_readout_dz) / 2.0 - clamp_plate_width),
406  single_tower_replica_logic,
407  "replicated_layers_placed",
408  single_tower_logic,
409  0, 0, OverlapCheck());
410 
411  // place array of fibers inside absorber
412  if (nFibers > 0 && nFibers == 5 && fiber_diam > 0)
413  {
414  std::string fiberName = "single_fiber_scintillator_solid" + std::to_string(type);
415  G4VSolid* single_scintillator_solid = new G4Tubs(fiberName,
416  0.0, fiber_diam / 2.0, (m_TowerDz[2] + fiber_extra_length) / 2, 0.0, 2 * M_PI);
417 
418  /* create logical volumes for scintillator and absorber plates to place inside single_tower */
419  G4Material* material_WLSFiber = GetWLSFiberFEMCMaterial();
420 
421  std::string fiberLogicName = "hEcal_scintillator_fiber_logic" + std::to_string(type);
422  G4LogicalVolume* single_scintillator_logic = new G4LogicalVolume(single_scintillator_solid,
423  material_WLSFiber,
424  fiberLogicName,
425  0, 0, 0);
426  m_ScintiLogicalVolSet.insert(single_scintillator_logic);
427  GetDisplayAction()->AddVolume(single_scintillator_logic, "Fiber");
428 
429  std::string name_scintillator = m_TowerLogicNamePrefix + "_single_fiber_scintillator" + std::to_string(type);
430 
431  new G4PVPlacement(0, G4ThreeVector(0, 0, -fiber_extra_length / 2 + tower_readout_dz / 2 - clamp_plate_width),
432  single_scintillator_logic,
433  name_scintillator + "_center",
434  single_tower_logic,
435  0, 0, OverlapCheck());
436 
437  new G4PVPlacement(0, G4ThreeVector(-m_TowerDx[2] / 4.0, m_TowerDy[2] / 4.0, -fiber_extra_length / 2 + tower_readout_dz / 2 - clamp_plate_width),
438  single_scintillator_logic,
439  name_scintillator + "_tl",
440  single_tower_logic,
441  0, 0, OverlapCheck());
442 
443  new G4PVPlacement(0, G4ThreeVector(m_TowerDx[2] / 4.0, -m_TowerDy[2] / 4.0, -fiber_extra_length / 2 + tower_readout_dz / 2 - clamp_plate_width),
444  single_scintillator_logic,
445  name_scintillator + "_bl",
446  single_tower_logic,
447  0, 0, OverlapCheck());
448 
449  new G4PVPlacement(0, G4ThreeVector(m_TowerDx[2] / 4.0, m_TowerDy[2] / 4.0, -fiber_extra_length / 2 + tower_readout_dz / 2 - clamp_plate_width),
450  single_scintillator_logic,
451  name_scintillator + "_tr",
452  single_tower_logic,
453  0, 0, OverlapCheck());
454 
455  new G4PVPlacement(0, G4ThreeVector(-m_TowerDx[2] / 4.0, -m_TowerDy[2] / 4.0, -fiber_extra_length / 2 + tower_readout_dz / 2 - clamp_plate_width),
456  single_scintillator_logic,
457  name_scintillator + "_br",
458  single_tower_logic,
459  0, 0, OverlapCheck());
460  }
461 
462  if (Verbosity() > 0)
463  {
464  std::cout << "PHG4ForwardEcalDetector: Building logical volume for single tower done." << std::endl;
465  }
466 
467  return single_tower_logic;
468 }
469 
471 {
472  /* Loop over all tower positions in vector and place tower */
473  for (std::map<std::string, towerposition>::iterator iterator = m_TowerPositionMap.begin(); iterator != m_TowerPositionMap.end(); ++iterator)
474  {
475  if (Verbosity() > 0)
476  {
477  std::cout << "PHG4ForwardEcalDetector: Place tower " << iterator->first
478  << " idx_j = " << iterator->second.idx_j << ", idx_k = " << iterator->second.idx_k
479  << " at x = " << iterator->second.x << " , y = " << iterator->second.y << " , z = " << iterator->second.z << std::endl;
480  }
481 
482  assert(iterator->second.type >= 0 && iterator->second.type <= 6);
483  G4LogicalVolume* singletower = singletowerIn[iterator->second.type];
484  int copyno = (iterator->second.idx_j << 16) + iterator->second.idx_k;
485 
486  G4PVPlacement* tower_placement =
487  new G4PVPlacement(0, G4ThreeVector(iterator->second.x, iterator->second.y, iterator->second.z),
488  singletower,
489  iterator->first,
490  ecalenvelope,
491  0, copyno, OverlapCheck());
492 
493  m_GdmlConfig->exclude_physical_vol(tower_placement);
494  }
495 
496  return 0;
497 }
498 
499 //_______________________________________________________________________
501 {
502  G4double density;
503  G4int ncomponents;
504  G4Material* material_ScintFEMC = new G4Material("PolystyreneFEMC", density = 1.03 * g / cm3, ncomponents = 2);
505  material_ScintFEMC->AddElement(GetDetectorElement("C"), 8);
506  material_ScintFEMC->AddElement(GetDetectorElement("H"), 8);
507 
508  if (m_doLightProp)
509  {
510  const G4int ntab = 4;
511 
512  G4double wls_Energy[] = {2.00 * eV, 2.87 * eV, 2.90 * eV,
513  3.47 * eV};
514 
515  G4double rIndexPstyrene[] = {1.5, 1.5, 1.5, 1.5};
516  G4double absorption1[] = {2. * cm, 2. * cm, 2. * cm, 2. * cm};
517  G4double scintilFast[] = {0.0, 0.0, 1.0, 1.0};
519  fMPTPStyrene->AddProperty("RINDEX", wls_Energy, rIndexPstyrene, ntab);
520  fMPTPStyrene->AddProperty("ABSLENGTH", wls_Energy, absorption1, ntab);
521  fMPTPStyrene->AddProperty("SCINTILLATIONCOMPONENT1", wls_Energy, scintilFast, ntab);
522  fMPTPStyrene->AddConstProperty("SCINTILLATIONYIELD", 10. / keV);
523  fMPTPStyrene->AddConstProperty("RESOLUTIONSCALE", 1.0);
524  fMPTPStyrene->AddConstProperty("SCINTILLATIONTIMECONSTANT", 10. * ns);
525 
526  material_ScintFEMC->SetMaterialPropertiesTable(fMPTPStyrene);
527  }
528  // Set the Birks Constant for the Polystyrene scintillator
529  material_ScintFEMC->GetIonisation()->SetBirksConstant(0.126 * mm / MeV);
530 
531  return material_ScintFEMC;
532 }
533 
534 //_______________________________________________________________________
536 {
537  //--------------------------------------------------
538  // TiO2
539  //--------------------------------------------------
540  G4double density, fractionmass;
541  G4int ncomponents;
542  G4Material* material_TiO2 = new G4Material("TiO2_FEMC", density = 1.52 * g / cm3, ncomponents = 2);
543  material_TiO2->AddElement(GetDetectorElement("Ti"), 1);
544  material_TiO2->AddElement(GetDetectorElement("O"), 2);
545 
546  //--------------------------------------------------
547  // Scintillator Coating - 15% TiO2 and 85% polystyrene by weight.
548  //--------------------------------------------------
549  //Coating_FEMC (Glass + Epoxy)
550  density = 1.86 * g / cm3;
551  G4Material* Coating_FEMC = new G4Material("Coating_FEMC", density, ncomponents = 2);
552  Coating_FEMC->AddMaterial(GetDetectorMaterial("Epoxy"), fractionmass = 0.80);
553  Coating_FEMC->AddMaterial(material_TiO2, fractionmass = 0.20);
554 
555  return Coating_FEMC;
556 }
557 
558 //_____________________________________________________________________________
560 {
561  G4OpticalSurface* surface = new G4OpticalSurface("ScintWrapB1");
562 
563  new G4LogicalSkinSurface("CrystalSurfaceL", vol, surface);
564 
565  surface->SetType(dielectric_metal);
566  surface->SetFinish(polished);
567  surface->SetModel(glisur);
568 
569  //crystal optical surface
570 
571  //surface material
572  // const G4int ntab = 2;
573  // G4double opt_en[] = {1.551*eV, 3.545*eV}; // 350 - 800 nm
574  // G4double reflectivity[] = {0.8, 0.8};
575  // G4double efficiency[] = {0.9, 0.9};
577  // surfmat->AddProperty("REFLECTIVITY", opt_en, reflectivity, ntab);
578  // surfmat->AddProperty("EFFICIENCY", opt_en, efficiency, ntab);
579  surface->SetMaterialPropertiesTable(surfmat);
580  //csurf->DumpInfo();
581 
582 } //SurfaceTable
583 //_______________________________________________________________________
585 {
586  if (Verbosity() > 0)
587  {
588  std::cout << "PHG4ForwardEcalDetector: Making WLSFiberFEMC material..." << std::endl;
589  }
590 
591  G4double density;
592  G4int ncomponents;
593 
594  G4Material* material_WLSFiberFEMC = new G4Material("WLSFiberFEMC", density = 1.18 * g / cm3, ncomponents = 3);
595  material_WLSFiberFEMC->AddElement(GetDetectorElement("C"), 5);
596  material_WLSFiberFEMC->AddElement(GetDetectorElement("H"), 8);
597  material_WLSFiberFEMC->AddElement(GetDetectorElement("O"), 2);
598  if (m_doLightProp)
599  {
600  const G4int ntab = 4;
601  G4double wls_Energy[] = {2.00 * eV, 2.87 * eV, 2.90 * eV,
602  3.47 * eV};
603 
604  G4double RefractiveIndexFiber[] = {1.6, 1.6, 1.6, 1.6};
605  G4double AbsFiber[] = {9.0 * m, 9.0 * m, 0.1 * mm, 0.1 * mm};
606  G4double EmissionFib[] = {1.0, 1.0, 0.0, 0.0};
607  // Add entries into properties table
609  mptWLSfiber->AddProperty("RINDEX", wls_Energy, RefractiveIndexFiber, ntab);
610  mptWLSfiber->AddProperty("WLSABSLENGTH", wls_Energy, AbsFiber, ntab);
611  mptWLSfiber->AddProperty("WLSCOMPONENT", wls_Energy, EmissionFib, ntab);
612  mptWLSfiber->AddConstProperty("WLSTIMECONSTANT", 0.5 * ns);
613  material_WLSFiberFEMC->SetMaterialPropertiesTable(mptWLSfiber);
614  }
615  if (Verbosity() > 0)
616  {
617  std::cout << "PHG4ForwardEcalDetector: Making WLSFiberFEMC material done." << std::endl;
618  }
619 
620  return material_WLSFiberFEMC;
621 }
622 
624 {
625  /* Open the datafile, if it won't open return an error */
626  std::ifstream istream_mapping;
627  istream_mapping.open(m_Params->get_string_param("mapping_file"));
628  if (!istream_mapping.is_open())
629  {
630  std::cout << "ERROR in PHG4ForwardEcalDetector: Failed to open mapping file " << m_Params->get_string_param("mapping_file") << std::endl;
631  gSystem->Exit(1);
632  }
633 
634  /* loop over lines in file */
635  std::string line_mapping;
636  while (getline(istream_mapping, line_mapping))
637  {
638  /* Skip lines starting with / including a '#' */
639  if (line_mapping.find("#") != std::string::npos)
640  {
641  if (Verbosity() > 0)
642  {
643  std::cout << "PHG4ForwardEcalDetector: SKIPPING line in mapping file: " << line_mapping << std::endl;
644  }
645  continue;
646  }
647 
648  std::istringstream iss(line_mapping);
649 
650  /* If line starts with keyword Tower, add to tower positions */
651  if (line_mapping.find("Tower ") != std::string::npos)
652  {
653  unsigned idx_j, idx_k, idx_l;
654  G4double pos_x, pos_y, pos_z;
655  G4double size_x, size_y, size_z;
656  G4double rot_x, rot_y, rot_z;
657  int type;
658  std::string dummys;
659 
660  /* read string- break if error */
661  if (!(iss >> dummys >> type >> idx_j >> idx_k >> idx_l >> pos_x >> pos_y >> pos_z >> size_x >> size_y >> size_z >> rot_x >> rot_y >> rot_z))
662  {
663  std::cout << "ERROR in PHG4ForwardEcalDetector: Failed to read line in mapping file " << m_Params->get_string_param("mapping_file") << std::endl;
664  gSystem->Exit(1);
665  }
666 
667  /* Construct unique name for tower */
668  /* Mapping file uses cm, this class uses mm for length */
669  std::ostringstream towername;
670  towername << m_TowerLogicNamePrefix << "_t_" << type << "_j_" << idx_j << "_k_" << idx_k;
671  /* Add Geant4 units */
672  pos_x = pos_x * cm;
673  pos_y = pos_y * cm;
674  pos_z = pos_z * cm;
675 
676  /* insert tower into tower map */
677  towerposition tower_new;
678  tower_new.x = pos_x;
679  tower_new.y = pos_y;
680  tower_new.z = pos_z;
681  tower_new.idx_j = idx_j;
682  tower_new.idx_k = idx_k;
683  tower_new.type = type;
684  m_TowerPositionMap.insert(std::make_pair(towername.str(), tower_new));
685  }
686  else
687  {
688  /* If this line is not a comment and not a tower, save parameter as string / value. */
689  std::string parname;
690  double parval;
691 
692  /* read string- break if error */
693  if (!(iss >> parname >> parval))
694  {
695  std::cout << "ERROR in PHG4ForwardEcalDetector: Failed to read line in mapping file " << m_Params->get_string_param("mapping_file") << std::endl;
696  gSystem->Exit(1);
697  }
698 
699  m_GlobalParameterMap.insert(std::make_pair(parname, parval));
700  }
701  }
702  /* Update member variables for global parameters based on parsed parameter file */
703  std::map<std::string, double>::iterator parit;
704  std::ostringstream twr;
705  for (int i = 0; i < 7; i++)
706  {
707  twr.str("");
708  twr << "Gtower" << i << "_dx";
709  parit = m_GlobalParameterMap.find(twr.str());
710  m_TowerDx[i] = parit->second * cm;
711  twr.str("");
712  twr << "Gtower" << i << "_dy";
713  parit = m_GlobalParameterMap.find(twr.str());
714  m_TowerDy[i] = parit->second * cm;
715  twr.str("");
716  twr << "Gtower" << i << "_dz";
717  parit = m_GlobalParameterMap.find(twr.str());
718  m_TowerDz[i] = parit->second * cm;
719  }
720  std::ostringstream rad;
721  for (int i = 0; i < 2; i++)
722  {
723  int index = i + 1;
724  rad.str("");
725  rad << "Gr" << index << "_inner";
726  parit = m_GlobalParameterMap.find(rad.str());
727  if (parit != m_GlobalParameterMap.end())
728  {
729  m_RMin[i] = parit->second * cm;
730  }
731  rad.str("");
732  rad << "Gr" << index << "_outer";
733  parit = m_GlobalParameterMap.find(rad.str());
734  if (parit != m_GlobalParameterMap.end())
735  {
736  m_RMax[i] = parit->second * cm;
737  }
738  }
739  parit = m_GlobalParameterMap.find("Gdz");
740  if (parit != m_GlobalParameterMap.end())
741  {
742  m_dZ = parit->second * cm;
743  }
744  parit = m_GlobalParameterMap.find("Gx0");
745  if (parit != m_GlobalParameterMap.end())
746  {
747  m_PlaceX = parit->second * cm;
748  }
749  parit = m_GlobalParameterMap.find("Gy0");
750  if (parit != m_GlobalParameterMap.end())
751  {
752  m_PlaceY = parit->second * cm;
753  }
754  parit = m_GlobalParameterMap.find("Gz0");
755  if (parit != m_GlobalParameterMap.end())
756  {
757  m_PlaceZ = parit->second * cm;
758  }
759  parit = m_GlobalParameterMap.find("Grot_x");
760  if (parit != m_GlobalParameterMap.end())
761  {
762  m_XRot = parit->second;
763  }
764  parit = m_GlobalParameterMap.find("Grot_y");
765  if (parit != m_GlobalParameterMap.end())
766  {
767  m_YRot = parit->second;
768  }
769  parit = m_GlobalParameterMap.find("Grot_z");
770  if (parit != m_GlobalParameterMap.end())
771  {
772  m_ZRot = parit->second;
773  }
774  parit = m_GlobalParameterMap.find("tower_type");
775  if (parit != m_GlobalParameterMap.end())
776  {
777  m_TowerType = parit->second;
778  }
779 
780  parit = m_GlobalParameterMap.find("xoffset");
781  if (parit != m_GlobalParameterMap.end())
782  m_Params->set_double_param("xoffset", parit->second);
783 
784  parit = m_GlobalParameterMap.find("yoffset");
785  if (parit != m_GlobalParameterMap.end())
786  m_Params->set_double_param("yoffset", parit->second);
787 
788  parit = m_GlobalParameterMap.find("width_coating");
789  if (parit != m_GlobalParameterMap.end())
790  m_Params->set_double_param("width_coating", parit->second);
791 
792  parit = m_GlobalParameterMap.find("clamp_plate_width");
793  if (parit != m_GlobalParameterMap.end())
794  m_Params->set_double_param("clamp_plate_width", parit->second);
795 
796  parit = m_GlobalParameterMap.find("fiber_diam");
797  if (parit != m_GlobalParameterMap.end())
798  m_Params->set_double_param("fiber_diam", parit->second);
799 
800  parit = m_GlobalParameterMap.find("nFibers");
801  if (parit != m_GlobalParameterMap.end())
802  m_Params->set_int_param("nFibers", parit->second);
803 
804  return 0;
805 }
806 
807 void PHG4ForwardEcalDetector::SetTowerDimensions(double dx, double dy, double dz, int type)
808 {
809  assert(type >= 0 && type <= 6);
810  m_TowerDx[type] = dx;
811  m_TowerDy[type] = dy;
812  m_TowerDz[type] = dz;
813 }