EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4LFHcalDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4LFHcalDetector.cc
1 #include "PHG4LFHcalDetector.h"
2 #include "PHG4LFHcalDisplayAction.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/G4ExtrudedSolid.hh>
15 #include <Geant4/G4LogicalBorderSurface.hh>
16 #include <Geant4/G4LogicalSkinSurface.hh>
17 #include <Geant4/G4LogicalVolume.hh>
18 #include <Geant4/G4Material.hh>
19 #include <Geant4/G4OpticalSurface.hh>
20 #include <Geant4/G4PVPlacement.hh>
21 #include <Geant4/G4PVReplica.hh>
22 #include <Geant4/G4RotationMatrix.hh> // for G4RotationMatrix
23 #include <Geant4/G4SubtractionSolid.hh>
24 #include <Geant4/G4SystemOfUnits.hh>
25 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
26 #include <Geant4/G4Torus.hh>
27 #include <Geant4/G4Transform3D.hh> // for G4Transform3D
28 #include <Geant4/G4Trap.hh>
29 #include <Geant4/G4Tubs.hh>
30 #include <Geant4/G4TwoVector.hh>
31 #include <Geant4/G4Types.hh> // for G4double, G4int
32 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
33 
34 #include <TSystem.h>
35 
36 #include <cmath>
37 #include <cstdlib>
38 #include <fstream>
39 #include <iostream>
40 #include <sstream>
41 #include <utility> // for pair, make_pair
42 
43 class G4VSolid;
44 class PHCompositeNode;
45 
46 using namespace std;
47 
48 //_______________________________________________________________________
49 PHG4LFHcalDetector::PHG4LFHcalDetector(PHG4Subsystem* subsys, PHCompositeNode* Node, PHParameters* parameters, const std::string& dnam)
50  : PHG4Detector(subsys, Node, dnam)
51  , m_DisplayAction(dynamic_cast<PHG4LFHcalDisplayAction*>(subsys->GetDisplayAction()))
52  , m_Params(parameters)
53  , m_ActiveFlag(m_Params->get_int_param("active"))
54  , m_AbsorberActiveFlag(m_Params->get_int_param("absorberactive"))
55  , m_TowerLogicNamePrefix("hHcalTower")
56  , m_SuperDetector("NONE")
57  , m_doLightProp(false)
58 {
59  if (m_Params->get_string_param("mapping_file").empty())
60  {
61  cout << "ERROR in PHG4LFHcalDetector: No mapping file specified. Abort detector construction." << endl;
62  cout << "Please run set_string_param(\"mapping_file\", std::string filename ) first." << endl;
63  gSystem->Exit(1);
64  }
65  /* Read parameters for detector construction and mappign from file */
66  ParseParametersFromTable();
67 }
68 //_______________________________________________________________________
69 int PHG4LFHcalDetector::IsInLFHcal(G4VPhysicalVolume* volume) const
70 {
71  G4LogicalVolume* mylogvol = volume->GetLogicalVolume();
72  if (m_ActiveFlag)
73  {
74  if (m_ScintiLogicalVolSet.find(mylogvol) != m_ScintiLogicalVolSet.end())
75  {
76  return 1;
77  }
78  }
79 
80  if (m_AbsorberActiveFlag)
81  {
82  if (m_AbsorberLogicalVolSet.find(mylogvol) != m_AbsorberLogicalVolSet.end())
83  {
84  return -1;
85  }
86  }
87  return 0;
88 }
89 
90 //_______________________________________________________________________
91 void PHG4LFHcalDetector::ConstructMe(G4LogicalVolume* logicWorld)
92 {
93  if (Verbosity() > 0)
94  {
95  std::cout << "PHG4LFHcalDetector: Begin Construction" << std::endl;
96  }
97 
98  // /* Read parameters for detector construction and mappign from file */
99  // ParseParametersFromTable();
100 
101  /* Create the cone envelope = 'world volume' for the crystal calorimeter */
102  recoConsts* rc = recoConsts::instance();
103  G4Material* WorldMaterial = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
104 
105  G4VSolid* beampipe_cutout = new G4Cons("LFHCAL_beampipe_cutout",
106  0, m_Params->get_double_param("rMin1") * cm,
107  0, m_Params->get_double_param("rMin2") * cm,
108  m_Params->get_double_param("dz") * cm / 2.0,
109  0, 2 * M_PI);
110  G4VSolid* hcal_envelope_solid = new G4Cons("LFHCAL_envelope_solid_cutout",
111  0, m_Params->get_double_param("rMax1") * cm,
112  0, m_Params->get_double_param("rMax2") * cm,
113  m_Params->get_double_param("dz") * cm / 2.0,
114  0, 2 * M_PI);
115  hcal_envelope_solid = new G4SubtractionSolid(G4String("LFHCAL_envelope_solid"), hcal_envelope_solid, beampipe_cutout, 0, G4ThreeVector(m_Params->get_double_param("xoffset") * cm, m_Params->get_double_param("yoffset") * cm, 0.));
116 
117  G4LogicalVolume* hcal_envelope_log = new G4LogicalVolume(hcal_envelope_solid, WorldMaterial, "hLFHCAL_envelope", 0, 0, 0);
118 
119  m_DisplayAction->AddVolume(hcal_envelope_log, "LFHcalEnvelope");
120 
121  /* Define rotation attributes for envelope cone */
122  G4RotationMatrix hcal_rotm;
123  hcal_rotm.rotateX(m_Params->get_double_param("rot_x") * deg);
124  hcal_rotm.rotateY(m_Params->get_double_param("rot_y") * deg);
125  hcal_rotm.rotateZ(m_Params->get_double_param("rot_z") * deg);
126 
127  /* Place envelope cone in simulation */
128  string name_envelope = m_TowerLogicNamePrefix + "_envelope";
129 
130  new G4PVPlacement(G4Transform3D(hcal_rotm, G4ThreeVector(m_Params->get_double_param("place_x") * cm,
131  m_Params->get_double_param("place_y") * cm,
132  m_Params->get_double_param("place_z") * cm)),
133  hcal_envelope_log, name_envelope, logicWorld, 0, false, OverlapCheck());
134 
135  /* Construct single calorimeter tower */
136  G4LogicalVolume* singletower = ConstructTower();
137 
138  /* Place calorimeter tower within envelope */
139  PlaceTower(hcal_envelope_log, singletower);
140 
141  return;
142 }
143 
144 //_______________________________________________________________________
145 G4LogicalVolume*
146 PHG4LFHcalDetector::ConstructTower()
147 {
148  if (Verbosity() > 0)
149  {
150  std::cout << "PHG4LFHcalDetector: Build logical volume for single tower..." << std::endl;
151  }
152 
153  //**********************************************************************************************
154  /* read variables from external inputs */
155  //**********************************************************************************************
156  G4double TowerDx = m_Params->get_double_param("tower_dx") * cm;
157  G4double TowerDy = m_Params->get_double_param("tower_dy") * cm;
158  G4double TowerDz = m_Params->get_double_param("tower_dz") * cm;
159  G4double WlsDw = m_Params->get_double_param("wls_dw") * cm;
160  G4double thick_frame_width = m_Params->get_double_param("frame_width") * cm;
161  G4double thin_frame_width = 0;
162  if (thick_frame_width > 0)
163  thin_frame_width = 0.1 * mm;
164  // double width_coating = m_Params->get_double_param("width_coating") * cm;
165  G4double thickness_absorber = m_Params->get_double_param("thickness_absorber") * cm;
166  G4double thickness_scintillator = m_Params->get_double_param("thickness_scintillator") * cm;
167  G4int nlayers = TowerDz / (thickness_absorber + thickness_scintillator);
168  G4Material* material_scintillator = GetScintillatorMaterial(); //GetDetectorMaterial(m_Params->get_string_param("scintillator"));
169  G4Material* material_absorber = GetDetectorMaterial(m_Params->get_string_param("absorber"));
170  G4Material* material_absorber_W = GetDetectorMaterial(m_Params->get_string_param("absorber_W"));
171  G4Material* material_wls = GetWLSFiberMaterial(); //GetDetectorMaterial(m_Params->get_string_param("scintillator"));
172  int embed_fiber = m_Params->get_int_param("embed_fiber");
173  G4double fiber_thickness = 0.2 * mm;
174  // G4double fiber_thickness = 1.0*mm;
175  //**********************************************************************************************
176  /* create logical volume for single tower */
177  //**********************************************************************************************
178  recoConsts* rc = recoConsts::instance();
179  G4Material* WorldMaterial = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
180  G4VSolid* single_tower_solid = new G4Box("single_tower_solid",
181  TowerDx / 2.0,
182  TowerDy / 2.0,
183  TowerDz / 2.0);
184 
185  G4LogicalVolume* single_tower_logic = new G4LogicalVolume(single_tower_solid,
186  WorldMaterial,
187  "single_tower_logic",
188  0, 0, 0);
189 
190  G4double notch_length = 2.0 * cm;
191  std::vector<G4TwoVector> poligon = {
192  {(-TowerDx) / 2.0 + thin_frame_width, (-TowerDy) / 2.0 + thick_frame_width},
193  {(-TowerDx) / 2.0 + thin_frame_width, (TowerDy) / 2.0 - thick_frame_width},
194  {(TowerDx - WlsDw) / 2.0 - thin_frame_width, (TowerDy) / 2.0 - thick_frame_width},
195  {(TowerDx - WlsDw) / 2.0 - thin_frame_width, (-TowerDy + notch_length) / 2.0 + thick_frame_width},
196  {(TowerDx) / 2.0 - thin_frame_width, (-TowerDy + notch_length) / 2.0 + thick_frame_width},
197  {(TowerDx) / 2.0 - thin_frame_width, (-TowerDy) / 2.0 + thick_frame_width}};
198  //**********************************************************************************************
199  /* create logical and geometry volumes for minitower read-out unit */
200  //**********************************************************************************************
201  std::vector<G4ExtrudedSolid::ZSection> zsections_miniblock = {{-(thickness_absorber + thickness_scintillator) / 2.0, {0, 0}, 1.0}, {(thickness_absorber + thickness_scintillator) / 2.0, {0, 0}, 1.}};
202  G4VSolid* miniblock_solid = new G4ExtrudedSolid("miniblock_solid", poligon, zsections_miniblock);
203 
204  G4LogicalVolume* miniblock_logic = new G4LogicalVolume(miniblock_solid,
205  WorldMaterial,
206  "miniblock_logic",
207  0, 0, 0);
208  G4LogicalVolume* miniblock_W_logic = new G4LogicalVolume(miniblock_solid,
209  WorldMaterial,
210  "miniblock_W_logic",
211  0, 0, 0);
212  m_DisplayAction->AddVolume(miniblock_logic, "Invisible");
213  m_DisplayAction->AddVolume(miniblock_W_logic, "Invisible");
214  //**********************************************************************************************
215  /* create logical & geometry volumes for scintillator and absorber plates to place inside mini read-out unit */
216  //**********************************************************************************************
217  std::vector<G4ExtrudedSolid::ZSection> zsections_scintillator = {{-thickness_scintillator / 2, {0, 0}, 1.0}, {thickness_scintillator / 2, {0, 0}, 1.}};
218 
219  G4VSolid* solid_scintillator = new G4ExtrudedSolid("single_plate_scintillator", poligon, zsections_scintillator);
220 
221  G4VPhysicalVolume* physvol_fiber_loop_0 = nullptr;
222  G4VPhysicalVolume* physvol_fiber_loop_1 = nullptr;
223  G4VPhysicalVolume* physvol_fiber_straight_1 = nullptr;
224  G4LogicalVolume* logic_embed_fiber_loop = nullptr;
225  G4LogicalVolume* logic_embed_fiber_loop_2 = nullptr;
226  G4LogicalVolume* logic_embed_fiber_straight_1 = nullptr;
227  G4double cutout_margin = 0.1 * fiber_thickness;
228  G4double fiber_bending_R_1 = notch_length / 2.0 - (TowerDy / 2 - thick_frame_width - (0.9 * (TowerDx - WlsDw - 2 * thin_frame_width) / 2.0) - (2 * cutout_margin));
229  G4double fiber_bending_R = 1.0 * cm;
230  if (embed_fiber)
231  {
232  G4VSolid* solid_embed_fiber_loop = new G4Torus("solid_embed_fiber_loop",
233  0, fiber_thickness / 2.0,
234  (0.9 * (TowerDx - WlsDw - 2 * thin_frame_width) / 2.0) - (fiber_thickness + 2 * cutout_margin) / 2.0,
235  -0.4 * M_PI, 1.9 * M_PI);
236  G4VSolid* solid_embed_fiber_loop_2 = new G4Torus("solid_embed_fiber_loop_2",
237  0, fiber_thickness / 2.0,
238  (2 * fiber_bending_R_1) / 2.0 - (fiber_thickness + 2 * cutout_margin) / 2.0,
239  -0.5 * M_PI, 0.5 * M_PI);
240  G4VSolid* solid_embed_fiber_straight_1 = new G4Tubs("solid_embed_fiber_straight_1",
241  0,
242  fiber_thickness / 2.0,
243  // (TowerDx/2.0 - fiber_bending_R_1 - thin_frame_width) / 2.0,
244  ((TowerDx - WlsDw + fiber_thickness) / 2.0 - (fiber_bending_R_1)) / 2.0,
245  0, 2.0 * M_PI);
246  solid_scintillator = new G4SubtractionSolid(G4String("single_plate_scintillator_cu"),
247  solid_scintillator, solid_embed_fiber_loop,
248  0, G4ThreeVector(0, 0, -thickness_scintillator / 2 + (fiber_thickness + cutout_margin) / 2.0));
249  solid_scintillator = new G4SubtractionSolid(G4String("single_plate_scintillator_cu2"),
250  solid_scintillator, solid_embed_fiber_loop_2,
251  0, G4ThreeVector((TowerDx - WlsDw + fiber_thickness / 2.0) / 2.0 - (fiber_bending_R_1) + (thin_frame_width) + cutout_margin, -0.9 * (TowerDx - WlsDw - thin_frame_width) / 2.0 + (fiber_bending_R_1), -thickness_scintillator / 2 + (fiber_thickness + cutout_margin) / 2.0));
252  G4RotationMatrix* wls_rotm_fibrcu = new G4RotationMatrix();
253  wls_rotm_fibrcu->rotateY(90 * deg);
254  solid_scintillator = new G4SubtractionSolid(G4String("single_plate_scintillator_cu3"),
255  solid_scintillator, solid_embed_fiber_straight_1,
256  wls_rotm_fibrcu, G4ThreeVector(((TowerDx - WlsDw + fiber_thickness) / 2.0 - (fiber_bending_R_1)) / 2.0, //(TowerDx/2-((TowerDx - WlsDw+fiber_thickness/2.0)/2.0-(fiber_bending_R_1)+(thin_frame_width)+cutout_margin))/2+fiber_bending_R_1/2+(fiber_thickness),
257  -0.9 * (TowerDx - WlsDw - thin_frame_width) / 2.0 + (fiber_thickness + 2 * cutout_margin) / 2.0, -thickness_scintillator / 2 + (fiber_thickness + cutout_margin) / 2.0));
258  logic_embed_fiber_loop = new G4LogicalVolume(solid_embed_fiber_loop,
259  material_wls, "logic_embed_fiber_loop",
260  0, 0, 0);
261  logic_embed_fiber_loop_2 = new G4LogicalVolume(solid_embed_fiber_loop_2,
262  material_wls, "logic_embed_fiber_loop_2",
263  0, 0, 0);
264  logic_embed_fiber_straight_1 = new G4LogicalVolume(solid_embed_fiber_straight_1,
265  material_wls, "logic_embed_fiber_straight_1",
266  0, 0, 0);
267  // if(m_doLightProp){
268  // SurfaceTable(logic_embed_fiber_loop);
269  // SurfaceTable(logic_embed_fiber_loop_2);
270  // SurfaceTable(logic_embed_fiber_straight_1);
271  // }
272  physvol_fiber_loop_0 = new G4PVPlacement(0, G4ThreeVector(0, 0, (thickness_absorber) / 2. - thickness_scintillator / 2 + (fiber_thickness + cutout_margin) / 2.0),
273  logic_embed_fiber_loop,
274  "embed_fiber_loop_placed",
275  miniblock_logic,
276  0, 0, OverlapCheck());
277  physvol_fiber_loop_1 = new G4PVPlacement(0, G4ThreeVector(
278  // (TowerDx - WlsDw+fiber_thickness/2.0)/2.0-(fiber_bending_R_1),
279  (TowerDx - WlsDw + fiber_thickness / 2.0) / 2.0 - (fiber_bending_R_1) + (thin_frame_width) + cutout_margin,
280  // (TowerDx)/2.0-(fiber_bending_R_1) - thin_frame_width,
281  -0.9 * (TowerDx - WlsDw - thin_frame_width) / 2.0 + (fiber_bending_R_1), (thickness_absorber) / 2. - thickness_scintillator / 2 + (fiber_thickness + cutout_margin) / 2.0),
282  logic_embed_fiber_loop_2,
283  "embed_fiber_loop_2_placed",
284  miniblock_logic,
285  0, 0, OverlapCheck());
286  G4RotationMatrix* wls_rotm_fibr = new G4RotationMatrix();
287  wls_rotm_fibr->rotateY(90 * deg);
288  physvol_fiber_straight_1 = new G4PVPlacement(wls_rotm_fibr, G4ThreeVector(
289  // (TowerDx/2.0 - fiber_bending_R_1 - thin_frame_width) / 2.0,
290  ((TowerDx - WlsDw + fiber_thickness) / 2.0 - (fiber_bending_R_1)) / 2.0,
291  // (TowerDx/2-((TowerDx - WlsDw+fiber_thickness/2.0)/2.0-(fiber_bending_R_1)+(thin_frame_width)+cutout_margin))/2,
292  -0.9 * (TowerDx - WlsDw - thin_frame_width) / 2.0 + (fiber_thickness + 2 * cutout_margin) / 2.0, (thickness_absorber) / 2. - thickness_scintillator / 2 + (fiber_thickness + cutout_margin) / 2.0),
293  logic_embed_fiber_straight_1,
294  "embed_fiber_straight_1_placed",
295  miniblock_logic,
296  0, 0, OverlapCheck());
297  if (m_doLightProp)
298  {
299  MakeBoundary(physvol_fiber_loop_0, physvol_fiber_loop_1, true);
300  MakeBoundary(physvol_fiber_loop_1, physvol_fiber_straight_1, true);
301  }
302  m_DisplayAction->AddVolume(logic_embed_fiber_loop, "WLSfiber");
303  m_DisplayAction->AddVolume(logic_embed_fiber_loop_2, "WLSfiber");
304  m_DisplayAction->AddVolume(logic_embed_fiber_straight_1, "WLSfiber");
305  }
306  // G4VSolid* solid_absorber = new G4Box("single_plate_absorber_solid",
307  // (TowerDx - WlsDw) / 2.0,
308  // (TowerDy ) / 2.0,
309  // thickness_absorber / 2.0);
310 
311  std::vector<G4ExtrudedSolid::ZSection> zsections_absorber = {{-thickness_absorber / 2, {0, 0}, 1.0}, {thickness_absorber / 2, {0, 0}, 1.}};
312  G4VSolid* solid_absorber = new G4ExtrudedSolid("single_plate_absorber_solid", poligon, zsections_absorber);
313  G4LogicalVolume* logic_absorber = new G4LogicalVolume(solid_absorber,
314  material_absorber,
315  "single_plate_absorber_logic",
316  0, 0, 0);
317  m_AbsorberLogicalVolSet.insert(logic_absorber);
318  G4LogicalVolume* logic_absorber_W = new G4LogicalVolume(solid_absorber,
319  material_absorber_W,
320  "single_plate_absorber_W_logic",
321  0, 0, 0);
322  m_AbsorberLogicalVolSet.insert(logic_absorber_W);
323  G4LogicalVolume* logic_scint = new G4LogicalVolume(solid_scintillator,
324  material_scintillator,
325  "hLFHCAL_scintillator_plate_logic",
326  0, 0, 0);
327  m_ScintiLogicalVolSet.insert(logic_scint);
328  // if(m_doLightProp){
329  // SurfaceTable(logic_scint);
330  // }
331  m_DisplayAction->AddVolume(logic_absorber, "Absorber");
332  m_DisplayAction->AddVolume(logic_absorber_W, "Absorber_W");
333  m_DisplayAction->AddVolume(logic_scint, "Scintillator");
334  string name_absorber = m_TowerLogicNamePrefix + "_single_plate_absorber";
335  string name_absorber_W = m_TowerLogicNamePrefix + "_single_plate_absorber_W";
336  string name_scintillator = m_TowerLogicNamePrefix + "_single_plate_scintillator";
337  // place Steel absorber and scintillator in miniblock
338  new G4PVPlacement(0, G4ThreeVector(0, 0, -thickness_scintillator / 2),
339  logic_absorber,
340  name_absorber,
341  miniblock_logic,
342  0, 0, OverlapCheck());
343 
344  G4VPhysicalVolume* scintillator_placed = new G4PVPlacement(0, G4ThreeVector(0, 0, (thickness_absorber) / 2.),
345  logic_scint,
346  name_scintillator,
347  miniblock_logic,
348  0, 0, OverlapCheck());
349  if (m_doLightProp && embed_fiber)
350  {
351  MakeBoundary_Fiber_Scint(physvol_fiber_loop_0, scintillator_placed);
352  MakeBoundary_Fiber_Scint(physvol_fiber_loop_1, scintillator_placed);
353  MakeBoundary_Fiber_Scint(physvol_fiber_straight_1, scintillator_placed);
354  }
355  // place Tungsten absorber and scintillator in a separate miniblock
356  new G4PVPlacement(0, G4ThreeVector(0, 0, -thickness_scintillator / 2),
357  logic_absorber_W,
358  name_absorber,
359  miniblock_W_logic,
360  0, 0, OverlapCheck());
361 
362  new G4PVPlacement(0, G4ThreeVector(0, 0, (thickness_absorber) / 2.),
363  logic_scint,
364  name_scintillator,
365  miniblock_W_logic,
366  0, 0, OverlapCheck());
367  if (embed_fiber)
368  {
369  G4VPhysicalVolume* physvol_fiber_loop_W_0 = new G4PVPlacement(0, G4ThreeVector(0, 0, (thickness_absorber) / 2. - thickness_scintillator / 2 + (fiber_thickness + cutout_margin) / 2.0),
370  logic_embed_fiber_loop,
371  "embed_fiber_loop_W_placed",
372  miniblock_W_logic,
373  0, 0, OverlapCheck());
374  G4VPhysicalVolume* physvol_fiber_loop_W_1 = new G4PVPlacement(0, G4ThreeVector((TowerDx - WlsDw + fiber_thickness / 2.0) / 2.0 - (fiber_bending_R_1) + (thin_frame_width) + cutout_margin,
375  // (TowerDx)/2.0-(fiber_bending_R_1) - thin_frame_width,
376  -0.9 * (TowerDx - WlsDw - thin_frame_width) / 2.0 + (fiber_bending_R_1), (thickness_absorber) / 2. - thickness_scintillator / 2 + (fiber_thickness + cutout_margin) / 2.0),
377  logic_embed_fiber_loop_2,
378  "embed_fiber_loop_W_2_placed",
379  miniblock_W_logic,
380  0, 0, OverlapCheck());
381  G4RotationMatrix* wls_rotm_fibr = new G4RotationMatrix();
382  wls_rotm_fibr->rotateY(90 * deg);
383  G4VPhysicalVolume* physvol_fiber_loop_W_2 = new G4PVPlacement(wls_rotm_fibr, G4ThreeVector(((TowerDx - WlsDw + fiber_thickness) / 2.0 - (fiber_bending_R_1)) / 2.0,
384  // (TowerDx/2-((TowerDx - WlsDw+fiber_thickness/2.0)/2.0-(fiber_bending_R_1)+(thin_frame_width)+cutout_margin))/2,
385  -0.9 * (TowerDx - WlsDw - thin_frame_width) / 2.0 + (fiber_thickness + 2 * cutout_margin) / 2.0, (thickness_absorber) / 2. - thickness_scintillator / 2 + (fiber_thickness + cutout_margin) / 2.0),
386  logic_embed_fiber_straight_1,
387  "embed_fiber_straight_W_1_placed",
388  miniblock_W_logic,
389  0, 0, OverlapCheck());
390  if (m_doLightProp)
391  {
392  MakeBoundary(physvol_fiber_loop_W_0, physvol_fiber_loop_W_1, true);
393  MakeBoundary(physvol_fiber_loop_W_1, physvol_fiber_loop_W_2, true);
394  }
395  }
396  //**********************************************************************************************
397  /* create wavelength shifting "fiber-block" */
398  //**********************************************************************************************
399  if (thin_frame_width > 0)
400  {
401  G4Material* material_frame = GetDetectorMaterial("G4_Fe");
402 
403  G4VSolid* solid_frame_plate = new G4Box("single_plate_frame",
404  (thin_frame_width) / 2, (TowerDy - 2 * thick_frame_width) / 2, TowerDz / 2);
405  G4LogicalVolume* logic_frame = new G4LogicalVolume(solid_frame_plate,
406  material_frame,
407  "hLFHCAL_frame_plate_logic",
408  0, 0, 0);
409  m_DisplayAction->AddVolume(logic_frame, "Frame");
410 
411  new G4PVPlacement(0, G4ThreeVector(-TowerDx / 2 + thin_frame_width / 2, 0, 0.),
412  logic_frame,
413  m_TowerLogicNamePrefix + "_single_plate_frame_left",
414  single_tower_logic,
415  0, 0, OverlapCheck());
416  new G4PVPlacement(0, G4ThreeVector(TowerDx / 2 - thin_frame_width / 2, 0, 0.),
417  logic_frame,
418  m_TowerLogicNamePrefix + "_single_plate_frame_right",
419  single_tower_logic,
420  0, 0, OverlapCheck());
421 
422  G4VSolid* solid_cover_plate = new G4Box("single_plate_cover",
423  TowerDx / 2, thick_frame_width / 2, TowerDz / 2);
424  G4LogicalVolume* logic_cover = new G4LogicalVolume(solid_cover_plate,
425  material_frame,
426  "hLFHCAL_cover_plate_logic",
427  0, 0, 0);
428  m_DisplayAction->AddVolume(logic_cover, "Frame");
429 
430  new G4PVPlacement(0, G4ThreeVector(0, -TowerDy / 2 + thick_frame_width / 2, 0.),
431  logic_cover,
432  m_TowerLogicNamePrefix + "_single_plate_cover_top",
433  single_tower_logic,
434  0, 0, OverlapCheck());
435  new G4PVPlacement(0, G4ThreeVector(0, TowerDy / 2 - thick_frame_width / 2, 0.),
436  logic_cover,
437  m_TowerLogicNamePrefix + "_single_plate_cover_bottom",
438  single_tower_logic,
439  0, 0, OverlapCheck());
440  }
441 
442  //**********************************************************************************************
443  /* create logical volume for single tower */
444  //**********************************************************************************************
445  double SteelTowerLength = TowerDz;
446  double WTowerLength = 0.;
447  int nLayersSteel = nlayers;
448  int nLayersTungsten = 0;
449  if (m_Params->get_int_param("usetailcatcher"))
450  {
451  SteelTowerLength -= 10 * (thickness_absorber + thickness_scintillator);
452  nLayersSteel -= 10;
453  nLayersTungsten = 10;
454  WTowerLength = 10 * (thickness_absorber + thickness_scintillator);
455  std::cout << "using 10 layer tungsten tailcatcher in LFHCAL" << std::endl;
456  }
457 
458  std::vector<G4ExtrudedSolid::ZSection> zsections_steeltower = {{-(SteelTowerLength) / 2.0, {0, 0}, 1.0}, {(SteelTowerLength) / 2.0, {0, 0}, 1.}};
459  G4VSolid* single_tower_solidRep = new G4ExtrudedSolid("single_tower_solidRep", poligon, zsections_steeltower);
460 
461  G4LogicalVolume* single_tower_logicRep = new G4LogicalVolume(single_tower_solidRep,
462  WorldMaterial,
463  "single_tower_logicRep",
464  0, 0, 0);
465  string name_tower = m_TowerLogicNamePrefix;
466  new G4PVReplica(name_tower, miniblock_logic, single_tower_logicRep,
467  kZAxis, nLayersSteel, thickness_absorber + thickness_scintillator, 0);
468 
469  new G4PVPlacement(0, G4ThreeVector(0, 0, -WTowerLength / 2),
470  single_tower_logicRep,
471  name_tower,
472  single_tower_logic,
473  0, 0, OverlapCheck());
474 
475  m_DisplayAction->AddVolume(single_tower_logicRep, "SingleTower");
476  m_DisplayAction->AddVolume(single_tower_logic, "SingleTower");
477 
478  if (m_Params->get_int_param("usetailcatcher"))
479  {
480  std::vector<G4ExtrudedSolid::ZSection> zsections_wtower = {{-(WTowerLength) / 2.0, {0, 0}, 1.0}, {(WTowerLength) / 2.0, {0, 0}, 1.}};
481  G4VSolid* single_tower_W_solidRep = new G4ExtrudedSolid("single_tower_W_solidRep", poligon, zsections_wtower);
482 
483  G4LogicalVolume* single_tower_W_logicRep = new G4LogicalVolume(single_tower_W_solidRep,
484  WorldMaterial,
485  "single_tower_W_logicRep",
486  0, 0, 0);
487  string name_tower_W = m_TowerLogicNamePrefix + "_W";
488  new G4PVReplica(name_tower_W, miniblock_W_logic, single_tower_W_logicRep,
489  kZAxis, nLayersTungsten, thickness_absorber + thickness_scintillator, 0);
490 
491  new G4PVPlacement(0, G4ThreeVector(0, 0, SteelTowerLength / 2),
492  single_tower_W_logicRep,
493  name_tower_W,
494  single_tower_logic,
495  0, 0, OverlapCheck());
496 
497  m_DisplayAction->AddVolume(single_tower_W_logicRep, "SingleTower_W");
498  }
499  if (embed_fiber)
500  {
501  G4double spacer_width = 1.5 * mm;
502  G4double add_spacing = spacer_width + 0.1 * mm;
503  G4double extraspacing = add_spacing;
504  for (int ilay = 0; ilay < nlayers; ilay++)
505  {
506  if ((SteelTowerLength + WTowerLength - ilay * (thickness_absorber + thickness_scintillator) - thickness_absorber - thickness_scintillator / 2 - fiber_bending_R / 2) / 2.0 > 0)
507  {
508  G4VSolid* solid_long_fiber_tmp = new G4Tubs("solid_long_fiber_tmp_" + std::to_string(ilay),
509  0,
510  fiber_thickness / 2.0,
511  (SteelTowerLength + WTowerLength - ilay * (thickness_absorber + thickness_scintillator) - thickness_absorber - fiber_thickness / 2 - fiber_bending_R / 2) / 2.0,
512  0, 2.0 * M_PI);
513 
514  G4LogicalVolume* logic_long_fiber_tmp = new G4LogicalVolume(solid_long_fiber_tmp,
515  material_wls,
516  "logic_long_fiber_tmp" + std::to_string(ilay),
517  0, 0, 0);
518  new G4PVPlacement(0, G4ThreeVector((TowerDx - WlsDw + fiber_thickness / 2.0) / 2.0, TowerDy / 2 - thick_frame_width - fiber_thickness - ilay * (1.01 * fiber_thickness) - extraspacing, (ilay * (thickness_absorber + thickness_scintillator) + thickness_absorber + fiber_thickness / 2 + fiber_bending_R / 2) / 2.0),
519  logic_long_fiber_tmp,
520  "phys_long_fiber_tmp" + std::to_string(ilay),
521  single_tower_logic,
522  0, 0, OverlapCheck());
523  m_DisplayAction->AddVolume(logic_long_fiber_tmp, "WLSfiber");
524  }
525 
526  // short fiber piece from loop to backward readout
527  G4double shortup_fiber_length = TowerDy / 2 - thick_frame_width - fiber_thickness - ilay * (1.01 * fiber_thickness) - extraspacing - fiber_bending_R / 2 + TowerDy / 2 - thick_frame_width - notch_length / 2;
528  G4VSolid* solid_shortup_fiber_tmp = new G4Tubs("solid_shortup_fiber_tmp_" + std::to_string(ilay),
529  0,
530  fiber_thickness / 2.0,
531  shortup_fiber_length / 2.0,
532  0, 2.0 * M_PI);
533  G4LogicalVolume* logic_shortup_fiber_tmp = new G4LogicalVolume(solid_shortup_fiber_tmp,
534  material_wls,
535  "logic_shortup_fiber_tmp" + std::to_string(ilay),
536  0, 0, 0);
537 
538  G4RotationMatrix* wls_rotm_fibr3 = new G4RotationMatrix();
539  wls_rotm_fibr3->rotateX(90 * deg);
540  new G4PVPlacement(wls_rotm_fibr3, G4ThreeVector((TowerDx - WlsDw + fiber_thickness / 2.0) / 2.0, TowerDy / 2 - thick_frame_width - fiber_thickness - ilay * (1.01 * fiber_thickness) - extraspacing - fiber_bending_R / 2 - (shortup_fiber_length / 2.0), -TowerDz / 2 + ilay * (thickness_absorber + thickness_scintillator) + thickness_absorber + fiber_thickness / 2),
541  logic_shortup_fiber_tmp,
542  "phys_shortup_fiber_tmp" + std::to_string(ilay),
543  single_tower_logic,
544  0, 0, OverlapCheck());
545  m_DisplayAction->AddVolume(logic_shortup_fiber_tmp, "WLSfiber");
546 
547  // curved fiber piece to backward readout
548  G4VSolid* solid_fiber_curved = new G4Torus("solid_fiber_curved",
549  0, fiber_thickness / 2.0,
550  (fiber_bending_R) / 2.0,
551  0.5 * M_PI, 0.5 * M_PI);
552  if (ilay == nlayers - 1)
553  {
554  solid_fiber_curved = new G4Torus("solid_fiber_curved_lastlayer",
555  0, fiber_thickness / 2.0,
556  (fiber_bending_R) / 2.0,
557  0.6 * M_PI, 0.4 * M_PI);
558  }
559  G4LogicalVolume* logic_fiber_curved = new G4LogicalVolume(solid_fiber_curved,
560  material_wls,
561  "logic_fiber_curved" + std::to_string(ilay),
562  0, 0, 0);
563 
564  G4RotationMatrix* wls_rotm_fibr2 = new G4RotationMatrix();
565  wls_rotm_fibr2->rotateY(90 * deg);
566  new G4PVPlacement(wls_rotm_fibr2, G4ThreeVector((TowerDx - WlsDw + fiber_thickness / 2.0) / 2.0, TowerDy / 2 - thick_frame_width - fiber_thickness - ilay * (1.01 * fiber_thickness) - extraspacing - fiber_bending_R / 2, -TowerDz / 2 + ilay * (thickness_absorber + thickness_scintillator) + thickness_absorber + fiber_thickness / 2 + fiber_bending_R / 2),
567  logic_fiber_curved,
568  "phys_fiber_curved" + std::to_string(ilay),
569  single_tower_logic,
570  0, 0, OverlapCheck());
571  m_DisplayAction->AddVolume(logic_fiber_curved, "WLSfiber");
572  if ((ilay + 1) % 10 == 0)
573  {
574  if (((SteelTowerLength + WTowerLength - ilay * (thickness_absorber + thickness_scintillator) - thickness_absorber - thickness_scintillator / 2 - thickness_absorber / 2) / 2.0) > 0)
575  {
576  G4VSolid* solid_spacer_tmp = new G4Box("solid_spacer_tmp_" + std::to_string(ilay),
577  (WlsDw) / 4,
578  spacer_width / 2.0,
579  (SteelTowerLength + WTowerLength - ilay * (thickness_absorber + thickness_scintillator) - thickness_absorber - thickness_scintillator / 2 - thickness_absorber / 2) / 2.0);
580 
581  G4LogicalVolume* logic_spacer_tmp = new G4LogicalVolume(solid_spacer_tmp,
582  GetDetectorMaterial("CFRP_INTT"), // carbon fiber + epoxy
583  "logic_spacer_tmp" + std::to_string(ilay),
584  0, 0, 0);
585  new G4PVPlacement(0, G4ThreeVector((TowerDx - WlsDw / 2) / 2.0 - thin_frame_width, TowerDy / 2 - thick_frame_width - fiber_thickness - (ilay + 0.5) * (1.01 * fiber_thickness) - extraspacing - add_spacing / 2, (ilay * (thickness_absorber + thickness_scintillator) + thickness_absorber + thickness_scintillator / 2 + thickness_absorber / 2) / 2.0),
586  logic_spacer_tmp,
587  "phys_spacer_tmp" + std::to_string(ilay),
588  single_tower_logic,
589  0, 0, OverlapCheck());
590  m_DisplayAction->AddVolume(logic_spacer_tmp, "Spacer");
591  }
592  extraspacing += add_spacing;
593  }
594  }
595  }
596  if (Verbosity() > 0)
597  {
598  std::cout << "PHG4LFHcalDetector: Building logical volume for single tower done." << std::endl;
599  }
600 
601  return single_tower_logic;
602 }
603 
604 int PHG4LFHcalDetector::PlaceTower(G4LogicalVolume* hcalenvelope, G4LogicalVolume* singletower)
605 {
606  /* Loop over all tower positions in vector and place tower */
607  for (std::map<std::string, towerposition>::iterator iterator = m_TowerPostionMap.begin(); iterator != m_TowerPostionMap.end(); ++iterator)
608  {
609  if (Verbosity() > 0)
610  {
611  std::cout << "PHG4LFHcalDetector: Place tower " << iterator->first
612  << " idx_j = " << iterator->second.idx_j << ", idx_k = " << iterator->second.idx_k
613  << " at x = " << iterator->second.x << " , y = " << iterator->second.y << " , z = " << iterator->second.z << std::endl;
614  }
615 
616  int copyno = (iterator->second.idx_j << 16) + iterator->second.idx_k;
617  // std::cout << "tower " << " idx_j = " << iterator->second.idx_j << ", idx_k = " << iterator->second.idx_k << " cpNo: " << copyno << std::endl;
618  new G4PVPlacement(0, G4ThreeVector(iterator->second.x, iterator->second.y, iterator->second.z),
619  singletower,
620  iterator->first,
621  hcalenvelope,
622  0, copyno, OverlapCheck());
623  }
624 
625  return 0;
626 }
627 
628 //_______________________________________________________________________
629 G4Material* PHG4LFHcalDetector::GetScintillatorMaterial()
630 {
631  G4double density;
632  G4int ncomponents;
633  G4Material* material_ScintFEMC = new G4Material("PolystyreneFEMC", density = 1.03 * g / cm3, ncomponents = 2);
634  material_ScintFEMC->AddElement(GetDetectorElement("C"), 8);
635  material_ScintFEMC->AddElement(GetDetectorElement("H"), 8);
636 
637  if (m_doLightProp)
638  {
639  // const G4int ntab = 31;
640  // G4double opt_en[] =
641  // { 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,
642  // 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,
643  // 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,
644  // 4.13281*eV }; // 350 - 800 nm
645  // G4double opt_abs[] =
646  // { 2.714*m, 3.619*m, 5.791*m, 4.343*m, 7.896*m, 5.429*m, 36.19*m, 17.37*m, 36.19*m, 5.429*m,
647  // 13.00*m, 14.50*m, 16.00*m, 18.00*m, 16.50*m, 17.00*m, 14.00*m, 16.00*m, 15.00*m, 14.50*m,
648  // 13.00*m, 12.00*m, 10.00*m, 8.000*m, 7.238*m, 4.000*m, 1.200*m, 0.500*m, 0.200*m, 0.200*m,
649  // 0.100*m };
650 
651  const G4int nEntries = 50;
652  G4double photonEnergy[nEntries] =
653  {2.00 * eV, 2.03 * eV, 2.06 * eV, 2.09 * eV, 2.12 * eV,
654  2.15 * eV, 2.18 * eV, 2.21 * eV, 2.24 * eV, 2.27 * eV,
655  2.30 * eV, 2.33 * eV, 2.36 * eV, 2.39 * eV, 2.42 * eV,
656  2.45 * eV, 2.48 * eV, 2.51 * eV, 2.54 * eV, 2.57 * eV,
657  2.60 * eV, 2.63 * eV, 2.66 * eV, 2.69 * eV, 2.72 * eV,
658  2.75 * eV, 2.78 * eV, 2.81 * eV, 2.84 * eV, 2.87 * eV,
659  2.90 * eV, 2.93 * eV, 2.96 * eV, 2.99 * eV, 3.02 * eV,
660  3.05 * eV, 3.08 * eV, 3.11 * eV, 3.14 * eV, 3.17 * eV,
661  3.20 * eV, 3.23 * eV, 3.26 * eV, 3.29 * eV, 3.32 * eV,
662  3.35 * eV, 3.38 * eV, 3.41 * eV, 3.44 * eV, 3.47 * eV};
663  G4double scintilFast[nEntries] =
664  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
665  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
666  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
667  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
668  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
669 
670  const G4int ntab = 4;
671 
672  G4double wls_Energy[] = {2.00 * eV, 2.87 * eV, 2.90 * eV,
673  3.47 * eV};
674 
675  G4double rIndexPstyrene[] = {1.5, 1.5, 1.5, 1.5};
676  G4double absorption1[] = {2. * cm, 2. * cm, 2. * cm, 2. * cm};
677  // G4double scintilFast[] = { 0.0, 0.0, 1.0, 1.0 };
678  G4MaterialPropertiesTable* fMPTPStyrene = new G4MaterialPropertiesTable();
679  fMPTPStyrene->AddProperty("RINDEX", wls_Energy, rIndexPstyrene, ntab);
680  fMPTPStyrene->AddProperty("ABSLENGTH", wls_Energy, absorption1, ntab);
681  // fMPTPStyrene->AddProperty("ABSLENGTH", opt_en, opt_abs, ntab);
682  // fMPTPStyrene->AddProperty("SCINTILLATIONCOMPONENT1", wls_Energy, scintilFast,ntab);
683  fMPTPStyrene->AddProperty("SCINTILLATIONCOMPONENT1", photonEnergy, scintilFast, nEntries);
684  fMPTPStyrene->AddConstProperty("SCINTILLATIONYIELD", 10. / keV);
685  fMPTPStyrene->AddConstProperty("RESOLUTIONSCALE", 1.0);
686  fMPTPStyrene->AddConstProperty("SCINTILLATIONTIMECONSTANT", 10. * ns);
687  // fMPTPStyrene->AddConstProperty("FASTTIMECONSTANT", 10.*ns);
688 
689  // fMPTPStyrene->AddProperty("FASTCOMPONENT",photonEnergy, scintilFast,nEntries);
690  material_ScintFEMC->SetMaterialPropertiesTable(fMPTPStyrene);
691  }
692  // Set the Birks Constant for the Polystyrene scintillator
693  material_ScintFEMC->GetIonisation()->SetBirksConstant(0.126 * mm / MeV);
694 
695  return material_ScintFEMC;
696 }
697 
698 //_______________________________________________________________________
699 G4Material* PHG4LFHcalDetector::GetCoatingMaterial()
700 {
701  //--------------------------------------------------
702  // TiO2
703  //--------------------------------------------------
704  G4double density, fractionmass;
705  G4int ncomponents;
706  G4Material* material_TiO2 = new G4Material("TiO2_FEMC", density = 1.52 * g / cm3, ncomponents = 2);
707  material_TiO2->AddElement(GetDetectorElement("Ti"), 1);
708  material_TiO2->AddElement(GetDetectorElement("O"), 2);
709 
710  //--------------------------------------------------
711  // Scintillator Coating - 15% TiO2 and 85% polystyrene by weight.
712  //--------------------------------------------------
713  //Coating_FEMC (Glass + Epoxy)
714  density = 1.86 * g / cm3;
715  G4Material* Coating_FEMC = new G4Material("Coating_FEMC", density, ncomponents = 2);
716  Coating_FEMC->AddMaterial(GetDetectorMaterial("Epoxy"), fractionmass = 0.80);
717  Coating_FEMC->AddMaterial(material_TiO2, fractionmass = 0.20);
718 
719  return Coating_FEMC;
720 }
721 
722 //_____________________________________________________________________________
723 void PHG4LFHcalDetector::SurfaceTable(G4LogicalVolume* vol)
724 {
725  G4OpticalSurface* surface = new G4OpticalSurface("ScintWrapB1");
726 
727  new G4LogicalSkinSurface("CrystalSurfaceL", vol, surface);
728 
729  surface->SetType(dielectric_metal);
730  surface->SetFinish(polished);
731  surface->SetModel(glisur);
732 
733  //crystal optical surface
734 
735  //surface material
736  const G4int ntab = 2;
737  G4double opt_en[] = {1.551 * eV, 3.545 * eV}; // 350 - 800 nm
738  G4double reflectivity[] = {0.9, 0.9};
739  G4double efficiency[] = {0.99, 0.99};
740  G4MaterialPropertiesTable* surfmat = new G4MaterialPropertiesTable();
741  surfmat->AddProperty("REFLECTIVITY", opt_en, reflectivity, ntab);
742  surfmat->AddProperty("EFFICIENCY", opt_en, efficiency, ntab);
743  surface->SetMaterialPropertiesTable(surfmat);
744  //csurf->DumpInfo();
745 
746 } //SurfaceTable
747 //_______________________________________________________________________
748 G4Material* PHG4LFHcalDetector::GetWLSFiberMaterial()
749 {
750  if (Verbosity() > 0)
751  {
752  std::cout << "PHG4ForwardEcalDetector: Making WLSFiberFEMC PMMA material..." << std::endl;
753  }
754 
755  G4double density;
756  G4int ncomponents;
757 
758  G4Material* material_WLSFiberFEMC = new G4Material("WLSFiberFEMC", density = 1.18 * g / cm3, ncomponents = 3);
759  material_WLSFiberFEMC->AddElement(GetDetectorElement("C"), 5);
760  material_WLSFiberFEMC->AddElement(GetDetectorElement("H"), 8);
761  material_WLSFiberFEMC->AddElement(GetDetectorElement("O"), 2);
762  if (m_doLightProp)
763  {
764  const G4int nEntries = 50;
765  G4double photonEnergy[nEntries] =
766  {2.00 * eV, 2.03 * eV, 2.06 * eV, 2.09 * eV, 2.12 * eV,
767  2.15 * eV, 2.18 * eV, 2.21 * eV, 2.24 * eV, 2.27 * eV,
768  2.30 * eV, 2.33 * eV, 2.36 * eV, 2.39 * eV, 2.42 * eV,
769  2.45 * eV, 2.48 * eV, 2.51 * eV, 2.54 * eV, 2.57 * eV,
770  2.60 * eV, 2.63 * eV, 2.66 * eV, 2.69 * eV, 2.72 * eV,
771  2.75 * eV, 2.78 * eV, 2.81 * eV, 2.84 * eV, 2.87 * eV,
772  2.90 * eV, 2.93 * eV, 2.96 * eV, 2.99 * eV, 3.02 * eV,
773  3.05 * eV, 3.08 * eV, 3.11 * eV, 3.14 * eV, 3.17 * eV,
774  3.20 * eV, 3.23 * eV, 3.26 * eV, 3.29 * eV, 3.32 * eV,
775  3.35 * eV, 3.38 * eV, 3.41 * eV, 3.44 * eV, 3.47 * eV};
776  G4double refractiveIndexWLSfiber[nEntries] =
777  {1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60,
778  1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60,
779  1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60,
780  1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60,
781  1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60};
782 
783  G4double absWLSfiber[nEntries] =
784  {5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m,
785  5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m,
786  5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m, 5.40 * m, 1.10 * m,
787  1.10 * m, 1.10 * m, 1.10 * m, 1.10 * m, 1.10 * m, 1.10 * m, 1. * mm, 1. * mm, 1. * mm, 1. * mm,
788  1. * mm, 1. * mm, 1. * mm, 1. * mm, 1. * mm, 1. * mm, 1. * mm, 1. * mm, 1. * mm, 1. * mm};
789 
790  G4double emissionFib[nEntries] =
791  {0.05, 0.10, 0.30, 0.50, 0.75, 1.00, 1.50, 1.85, 2.30, 2.75,
792  3.25, 3.80, 4.50, 5.20, 6.00, 7.00, 8.50, 9.50, 11.1, 12.4,
793  12.9, 13.0, 12.8, 12.3, 11.1, 11.0, 12.0, 11.0, 17.0, 16.9,
794  15.0, 9.00, 2.50, 1.00, 0.05, 0.00, 0.00, 0.00, 0.00, 0.00,
795  0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00};
796  // Add entries into properties table
797  G4MaterialPropertiesTable* mptWLSfiber = new G4MaterialPropertiesTable();
798  mptWLSfiber->AddProperty("RINDEX", photonEnergy, refractiveIndexWLSfiber, nEntries);
799  mptWLSfiber->AddProperty("WLSABSLENGTH", photonEnergy, absWLSfiber, nEntries);
800  mptWLSfiber->AddProperty("WLSCOMPONENT", photonEnergy, emissionFib, nEntries);
801  mptWLSfiber->AddConstProperty("WLSTIMECONSTANT", 0.5 * ns);
802  material_WLSFiberFEMC->SetMaterialPropertiesTable(mptWLSfiber);
803  }
804  if (Verbosity() > 0)
805  {
806  std::cout << "PHG4ForwardEcalDetector: Making WLSFiberFEMC material done." << std::endl;
807  }
808 
809  return material_WLSFiberFEMC;
810 }
811 
812 //_____________________________________________________________________________
813 void PHG4LFHcalDetector::MakeBoundary(G4VPhysicalVolume* crystal, G4VPhysicalVolume* opdet, bool isFiber)
814 {
815  //optical boundary between the crystal and optical photons detector
816 
817  G4OpticalSurface* surf = new G4OpticalSurface("OpDetS");
818  surf->SetType(dielectric_dielectric); // photons go to the detector, must have rindex defined
819  // surf->SetType(dielectric_metal); // photon is absorbed when reaching the detector, no material rindex required
820  //surf->SetFinish(ground);
821  surf->SetFinish(polished);
822  //surf->SetModel(unified);
823  surf->SetModel(glisur);
824 
825  new G4LogicalBorderSurface("OpDetB", crystal, opdet, surf);
826 
827  const G4int ntab = 2;
828  G4double opt_en[] = {1.551 * eV, 3.545 * eV}; // 350 - 800 nm
829  //G4double reflectivity[] = {0., 0.};
830  G4double reflectivityFiber[] = {0.0, 0.0};
831  //G4double reflectivity[] = {0.9, 0.9};
832  //G4double reflectivity[] = {1., 1.};
833  G4double efficiency[] = {1., 1.};
834  G4double RefractiveIndexBoundary[] = {1.0, 1.0};
835 
836  G4MaterialPropertiesTable* surfmat = new G4MaterialPropertiesTable();
837  surfmat->AddProperty("EFFICIENCY", opt_en, efficiency, ntab);
838  surfmat->AddProperty("RINDEX", opt_en, RefractiveIndexBoundary, ntab);
839  surfmat->AddProperty("REFLECTIVITY", opt_en, reflectivityFiber, ntab);
840  surf->SetMaterialPropertiesTable(surfmat);
841 
842 } //MakeBoundary
843 
844 //_____________________________________________________________________________
845 void PHG4LFHcalDetector::MakeBoundary_Fiber_Scint(G4VPhysicalVolume* fiber, G4VPhysicalVolume* scinti)
846 {
847  //optical boundary between the fiber and scintillator
848 
849  G4OpticalSurface* ScintToFiberS = new G4OpticalSurface("ScintToFiberS");
850  ScintToFiberS->SetType(dielectric_dielectric); // photons go to the detector, must have rindex defined
851  // ScintToFiberS->SetType(dielectric_metal); // photon is absorbed when reaching the detector, no material rindex required
852  //ScintToFiberS->SetFinish(ground);
853  ScintToFiberS->SetFinish(groundair);
854  // ScintToFiberS->SetFinish(polished);
855  //ScintToFiberS->SetModel(unified);
856  ScintToFiberS->SetModel(glisur);
857 
858  new G4LogicalBorderSurface("ScintToFiberB", scinti, fiber, ScintToFiberS);
859 
860  const G4int ntab = 2;
861  G4double opt_en[] = {1.551 * eV, 3.545 * eV}; // 350 - 800 nm
862  //G4double reflectivity[] = {0., 0.};
863  G4double reflectivity[] = {0.01, 0.01};
864  //G4double reflectivity[] = {0.9, 0.9};
865  //G4double reflectivity[] = {1., 1.};
866  G4double efficiency[] = {1., 1.};
867  G4double RefractiveIndexBoundary[] = {1.4, 1.4};
868 
869  G4MaterialPropertiesTable* ScintToFiberSmat = new G4MaterialPropertiesTable();
870  ScintToFiberSmat->AddProperty("EFFICIENCY", opt_en, efficiency, ntab);
871  ScintToFiberSmat->AddProperty("RINDEX", opt_en, RefractiveIndexBoundary, ntab);
872  ScintToFiberSmat->AddProperty("REFLECTIVITY", opt_en, reflectivity, ntab);
873  ScintToFiberS->SetMaterialPropertiesTable(ScintToFiberSmat);
874 
875  //optical boundary between the fiber and scintillator
876 
877  G4OpticalSurface* FiberToScintSurface = new G4OpticalSurface("ScintToFiberS");
878  // FiberToScintSurface->SetType(dielectric_dielectric); // photons go to the detector, must have rindex defined
879  FiberToScintSurface->SetType(dielectric_metal); // photon is absorbed when reaching the detector, no material rindex required
880  //FiberToScintSurface->SetFinish(ground);
881  // FiberToScintSurface->SetFinish(groundair);
882  FiberToScintSurface->SetFinish(polished);
883  //FiberToScintSurface->SetModel(unified);
884  FiberToScintSurface->SetModel(glisur);
885 
886  new G4LogicalBorderSurface("ScintToFiberB", scinti, fiber, FiberToScintSurface);
887 
888  const G4int ntab2 = 2;
889  G4double opt_en2[] = {1.551 * eV, 3.545 * eV}; // 350 - 800 nm
890  //G4double reflectivity[] = {0., 0.};
891  // G4double reflectivity2[] = {0.01, 0.01};
892  //G4double reflectivity[] = {0.9, 0.9};
893  G4double reflectivity2[] = {1., 1.};
894  G4double efficiency2[] = {1., 1.};
895  G4double RefractiveIndexBoundary2[] = {1.6, 1.6};
896 
897  G4MaterialPropertiesTable* FiberToScintSurfacemat = new G4MaterialPropertiesTable();
898  FiberToScintSurfacemat->AddProperty("EFFICIENCY", opt_en2, efficiency2, ntab2);
899  FiberToScintSurfacemat->AddProperty("RINDEX", opt_en2, RefractiveIndexBoundary2, ntab2);
900  FiberToScintSurfacemat->AddProperty("REFLECTIVITY", opt_en2, reflectivity2, ntab2);
901  FiberToScintSurface->SetMaterialPropertiesTable(FiberToScintSurfacemat);
902 
903 } //MakeBoundary
904 
905 int PHG4LFHcalDetector::ParseParametersFromTable()
906 {
907  /* Open the datafile, if it won't open return an error */
908  ifstream istream_mapping;
909  istream_mapping.open(m_Params->get_string_param("mapping_file"));
910  if (!istream_mapping.is_open())
911  {
912  std::cout << "ERROR in PHG4LFHcalDetector: Failed to open mapping file " << m_Params->get_string_param("mapping_file") << std::endl;
913  gSystem->Exit(1);
914  }
915 
916  /* loop over lines in file */
917  string line_mapping;
918  while (getline(istream_mapping, line_mapping))
919  {
920  /* Skip lines starting with / including a '#' */
921  if (line_mapping.find("#") != string::npos)
922  {
923  if (Verbosity() > 0)
924  {
925  std::cout << "PHG4LFHcalDetector: SKIPPING line in mapping file: " << line_mapping << std::endl;
926  }
927  continue;
928  }
929 
930  istringstream iss(line_mapping);
931 
932  /* If line starts with keyword Tower, add to tower positions */
933  if (line_mapping.find("Tower ") != string::npos)
934  {
935  unsigned idx_j, idx_k, idx_l;
936  G4double pos_x, pos_y, pos_z;
937  G4double size_x, size_y, size_z;
938  G4double rot_x, rot_y, rot_z;
939  G4double dummy;
940  string dummys;
941 
942  /* read string- break if error */
943  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))
944  {
945  cout << "ERROR in PHG4LFHcalDetector: Failed to read line in mapping file " << m_Params->get_string_param("mapping_file") << std::endl;
946  gSystem->Exit(1);
947  }
948 
949  /* Construct unique name for tower */
950  /* Mapping file uses cm, this class uses mm for length */
951  ostringstream towername;
952  towername.str("");
953  towername << m_TowerLogicNamePrefix << "_j_" << idx_j << "_k_" << idx_k;
954 
955  /* Add Geant4 units */
956  pos_x = pos_x * cm;
957  pos_y = pos_y * cm;
958  pos_z = pos_z * cm;
959 
960  /* insert tower into tower map */
961  towerposition tower_new;
962  tower_new.x = pos_x;
963  tower_new.y = pos_y;
964  tower_new.z = pos_z;
965  tower_new.idx_j = idx_j;
966  tower_new.idx_k = idx_k;
967  m_TowerPostionMap.insert(make_pair(towername.str(), tower_new));
968  }
969  else
970  {
971  /* If this line is not a comment and not a tower, save parameter as string / value. */
972  string parname;
973  G4double parval;
974 
975  /* read string- break if error */
976  if (!(iss >> parname >> parval))
977  {
978  cout << "ERROR in PHG4LFHcalDetector: Failed to read line in mapping file " << m_Params->get_string_param("mapping_file") << std::endl;
979  gSystem->Exit(1);
980  }
981 
982  m_GlobalParameterMap.insert(make_pair(parname, parval));
983  }
984  }
985 
986  /* Update member variables for global parameters based on parsed parameter file */
987  std::map<string, G4double>::iterator parit;
988 
989  parit = m_GlobalParameterMap.find("Gtower_dx");
990  if (parit != m_GlobalParameterMap.end())
991  {
992  m_Params->set_double_param("tower_dx", parit->second); // in cm
993  }
994 
995  parit = m_GlobalParameterMap.find("Gtower_dy");
996  if (parit != m_GlobalParameterMap.end())
997  {
998  m_Params->set_double_param("tower_dy", parit->second); // in cm
999  }
1000 
1001  parit = m_GlobalParameterMap.find("Gtower_dz");
1002  if (parit != m_GlobalParameterMap.end())
1003  {
1004  m_Params->set_double_param("tower_dz", parit->second); // in cm
1005  }
1006 
1007  parit = m_GlobalParameterMap.find("Gr1_inner");
1008  if (parit != m_GlobalParameterMap.end())
1009  {
1010  m_Params->set_double_param("rMin1", parit->second);
1011  }
1012 
1013  parit = m_GlobalParameterMap.find("Gr1_outer");
1014  if (parit != m_GlobalParameterMap.end())
1015  {
1016  m_Params->set_double_param("rMax1", parit->second);
1017  }
1018 
1019  parit = m_GlobalParameterMap.find("Gr2_inner");
1020  if (parit != m_GlobalParameterMap.end())
1021  {
1022  m_Params->set_double_param("rMin2", parit->second);
1023  }
1024 
1025  parit = m_GlobalParameterMap.find("Gr2_outer");
1026  if (parit != m_GlobalParameterMap.end())
1027  {
1028  m_Params->set_double_param("rMax2", parit->second);
1029  }
1030 
1031  parit = m_GlobalParameterMap.find("Gdz");
1032  if (parit != m_GlobalParameterMap.end())
1033  {
1034  m_Params->set_double_param("dz", parit->second);
1035  }
1036 
1037  parit = m_GlobalParameterMap.find("Gx0");
1038  if (parit != m_GlobalParameterMap.end())
1039  {
1040  m_Params->set_double_param("place_x", parit->second);
1041  }
1042 
1043  parit = m_GlobalParameterMap.find("Gy0");
1044  if (parit != m_GlobalParameterMap.end())
1045  {
1046  m_Params->set_double_param("place_y", parit->second);
1047  }
1048 
1049  parit = m_GlobalParameterMap.find("Gz0");
1050  if (parit != m_GlobalParameterMap.end())
1051  {
1052  m_Params->set_double_param("place_z", parit->second);
1053  }
1054 
1055  parit = m_GlobalParameterMap.find("Grot_x");
1056  if (parit != m_GlobalParameterMap.end())
1057  {
1058  m_Params->set_double_param("rot_x", parit->second * rad / deg);
1059  }
1060 
1061  parit = m_GlobalParameterMap.find("Grot_y");
1062  if (parit != m_GlobalParameterMap.end())
1063  {
1064  m_Params->set_double_param("rot_y", parit->second * rad / deg);
1065  }
1066 
1067  parit = m_GlobalParameterMap.find("Grot_z");
1068  if (parit != m_GlobalParameterMap.end())
1069  {
1070  m_Params->set_double_param("rot_z", parit->second * rad / deg);
1071  }
1072 
1073  parit = m_GlobalParameterMap.find("thickness_absorber");
1074  if (parit != m_GlobalParameterMap.end())
1075  m_Params->set_double_param("thickness_absorber", parit->second);
1076 
1077  parit = m_GlobalParameterMap.find("thickness_scintillator");
1078  if (parit != m_GlobalParameterMap.end())
1079  m_Params->set_double_param("thickness_scintillator", parit->second);
1080 
1081  parit = m_GlobalParameterMap.find("nlayerspertowerseg");
1082  if (parit != m_GlobalParameterMap.end())
1083  m_Params->set_int_param("nlayerspertowerseg", (int) parit->second);
1084 
1085  parit = m_GlobalParameterMap.find("xoffset");
1086  if (parit != m_GlobalParameterMap.end())
1087  m_Params->set_double_param("xoffset", parit->second);
1088 
1089  parit = m_GlobalParameterMap.find("yoffset");
1090  if (parit != m_GlobalParameterMap.end())
1091  m_Params->set_double_param("yoffset", parit->second);
1092 
1093  parit = m_GlobalParameterMap.find("frame_width");
1094  if (parit != m_GlobalParameterMap.end())
1095  m_Params->set_double_param("frame_width", parit->second);
1096 
1097  parit = m_GlobalParameterMap.find("width_coating");
1098  if (parit != m_GlobalParameterMap.end())
1099  m_Params->set_double_param("width_coating", parit->second);
1100 
1101  parit = m_GlobalParameterMap.find("usetailcatcher");
1102  if (parit != m_GlobalParameterMap.end())
1103  m_Params->set_int_param("usetailcatcher", parit->second);
1104 
1105  parit = m_GlobalParameterMap.find("embed_fiber");
1106  if (parit != m_GlobalParameterMap.end())
1107  m_Params->set_int_param("embed_fiber", parit->second);
1108 
1110  if (m_Params->get_int_param("usetailcatcher"))
1111  {
1112  m_Params->set_double_param("zdepthcatcheroffset", m_Params->get_double_param("place_z") + (m_Params->get_double_param("tower_dz") / 2) - (10 * (m_Params->get_double_param("thickness_scintillator") + m_Params->get_double_param("thickness_absorber"))));
1113  m_Params->set_int_param("nLayerOffsetTailcatcher", (int) ((m_Params->get_double_param("tower_dz") - (10 * (m_Params->get_double_param("thickness_scintillator") + m_Params->get_double_param("thickness_absorber")))) / (m_Params->get_double_param("thickness_scintillator") + m_Params->get_double_param("thickness_absorber"))));
1114  }
1115  if (Verbosity() > 1)
1116  {
1117  std::cout << "PHG4 detector LFHCal - Absorber: " << m_Params->get_double_param("thickness_absorber") << " cm\t Scintilator: " << m_Params->get_double_param("thickness_scintillator") << " cm\t layers per segment: " << m_Params->get_int_param("nlayerspertowerseg") << std::endl;
1118  }
1119 
1120  return 0;
1121 }