ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4BackwardHcalDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4BackwardHcalDetector.cc
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/G4LogicalVolume.hh>
15 #include <Geant4/G4Material.hh>
16 #include <Geant4/G4PVPlacement.hh>
17 #include <Geant4/G4RotationMatrix.hh> // for G4RotationMatrix
18 #include <Geant4/G4SystemOfUnits.hh>
19 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
20 #include <Geant4/G4Transform3D.hh> // for G4Transform3D
21 #include <Geant4/G4Types.hh> // for G4double, G4int
22 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
23 
24 #include <TSystem.h>
25 
26 #include <cmath>
27 #include <cstdlib>
28 #include <fstream>
29 #include <iostream>
30 #include <sstream>
31 #include <utility> // for pair, make_pair
32 
33 class G4VSolid;
34 class PHCompositeNode;
35 
36 using namespace std;
37 
38 //_______________________________________________________________________
40  : PHG4Detector(subsys, Node, dnam)
41  , m_DisplayAction(dynamic_cast<PHG4BackwardHcalDisplayAction*>(subsys->GetDisplayAction()))
42  , m_Params(parameters)
43  , m_ActiveFlag(m_Params->get_int_param("active"))
44  , m_AbsorberActiveFlag(m_Params->get_int_param("absorberactive"))
45  , m_TowerLogicNamePrefix("hHcalTower")
46  , m_SuperDetector("NONE")
47 {
48 }
49 //_______________________________________________________________________
51 {
52  G4LogicalVolume* mylogvol = volume->GetLogicalVolume();
53  if (m_ActiveFlag)
54  {
55  if (m_ScintiLogicalVolSet.find(mylogvol) != m_ScintiLogicalVolSet.end())
56  {
57  return 1;
58  }
59  }
60 
62  {
63  if (m_AbsorberLogicalVolSet.find(mylogvol) != m_AbsorberLogicalVolSet.end())
64  {
65  return -1;
66  }
67  }
68  return 0;
69 }
70 
71 //_______________________________________________________________________
73 {
74  if (Verbosity() > 0)
75  {
76  std::cout << "PHG4BackwardHcalDetector: Begin Construction" << std::endl;
77  }
78 
79  if (m_Params->get_string_param("mapping_file").empty())
80  {
81  cout << "ERROR in PHG4BackwardHcalDetector: No mapping file specified. Abort detector construction." << endl;
82  cout << "Please run set_string_param(\"mapping_file\", std::string filename ) first." << endl;
83  gSystem->Exit(1);
84  }
85 
86  /* Read parameters for detector construction and mappign from file */
88 
89  /* Create the cone envelope = 'world volume' for the crystal calorimeter */
92 
93  G4VSolid* hcal_envelope_solid = new G4Cons("hHcal_envelope_solid",
94  m_Params->get_double_param("rMin1") * cm,
95  m_Params->get_double_param("rMax1") * cm,
96  m_Params->get_double_param("rMin2") * cm,
97  m_Params->get_double_param("rMax2") * cm,
98  m_Params->get_double_param("dz") * cm / 2.,
99  0., 2. * M_PI);
100 
101  G4LogicalVolume* hcal_envelope_log = new G4LogicalVolume(hcal_envelope_solid, WorldMaterial, "hHcal_envelope", 0, 0, 0);
102 
103  m_DisplayAction->AddVolume(hcal_envelope_log, "FHcalEnvelope");
104 
105  /* Define rotation attributes for envelope cone */
106  G4RotationMatrix hcal_rotm;
107  hcal_rotm.rotateX(m_Params->get_double_param("rot_x") * deg);
108  hcal_rotm.rotateY(m_Params->get_double_param("rot_y") * deg);
109  hcal_rotm.rotateZ(m_Params->get_double_param("rot_z") * deg);
110 
111  /* Place envelope cone in simulation */
112  string name_envelope = m_TowerLogicNamePrefix + "_envelope";
113 
115  m_Params->get_double_param("place_y") * cm,
116  m_Params->get_double_param("place_z") * cm)),
117  hcal_envelope_log, name_envelope, logicWorld, 0, false, OverlapCheck());
118 
119  /* Construct single calorimeter tower */
120  G4LogicalVolume* singletower = ConstructTower();
121 
122  /* Place calorimeter tower within envelope */
123  PlaceTower(hcal_envelope_log, singletower);
124 
125  return;
126 }
127 
128 //_______________________________________________________________________
131 {
132  if (Verbosity() > 0)
133  {
134  std::cout << "PHG4BackwardHcalDetector: Build logical volume for single tower..." << std::endl;
135  }
136 
137  /* create logical volume for single tower */
140  double TowerDx = m_Params->get_double_param("tower_dx") * cm;
141  double TowerDy = m_Params->get_double_param("tower_dy") * cm;
142  double TowerDz = m_Params->get_double_param("tower_dz") * cm;
143  double WlsDw = m_Params->get_double_param("wls_dw") * cm;
144  double SupportDw = m_Params->get_double_param("support_dw") * cm;
145  G4VSolid* single_tower_solid = new G4Box("single_tower_solid",
146  TowerDx / 2.0,
147  TowerDy / 2.0,
148  TowerDz / 2.0);
149 
150  G4LogicalVolume* single_tower_logic = new G4LogicalVolume(single_tower_solid,
151  WorldMaterial,
152  "single_tower_logic",
153  0, 0, 0);
154 
155  /* create geometry volumes for scintillator and absorber plates to place inside single_tower */
156  // based on STAR forward upgrade design: https://drupal.star.bnl.gov/STAR/files/ForwardUpgrade.v20.pdf
157  G4double thickness_absorber = m_Params->get_double_param("thickness_absorber") * cm;
158  G4double thickness_scintillator = m_Params->get_double_param("thickness_scintillator") * cm;
159  G4int nlayers = TowerDz / (thickness_absorber + thickness_scintillator);
160 
161  G4VSolid* solid_absorber = new G4Box("single_plate_absorber_solid",
162  (TowerDx - WlsDw) / 2.0,
163  (TowerDy - SupportDw) / 2.0,
164  thickness_absorber / 2.0);
165 
166  G4VSolid* solid_scintillator = new G4Box("single_plate_scintillator",
167  (TowerDx - WlsDw) / 2.0,
168  (TowerDy - SupportDw) / 2.0,
169  thickness_scintillator / 2.0);
170 
171  G4VSolid* solid_WLS_plate = new G4Box("single_plate_wls",
172  (WlsDw) / 2.0,
173  (TowerDy - SupportDw) / 2.0,
174  TowerDz / 2.0);
175 
176  G4VSolid* solid_support_plate = new G4Box("single_plate_support",
177  (TowerDx) / 2.0,
178  (SupportDw) / 2.0,
179  TowerDz / 2.0);
180 
181  /* create logical volumes for scintillator and absorber plates to place inside single_tower */
182  G4Material* material_scintillator = GetDetectorMaterial(m_Params->get_string_param("scintillator"));
183  G4Material* material_absorber = GetDetectorMaterial(m_Params->get_string_param("absorber"));
184  G4Material* material_wls = GetDetectorMaterial(m_Params->get_string_param("scintillator"));
185  G4Material* material_support = GetDetectorMaterial(m_Params->get_string_param("support"));
186 
187  G4LogicalVolume* logic_absorber = new G4LogicalVolume(solid_absorber,
188  material_absorber,
189  "single_plate_absorber_logic",
190  0, 0, 0);
191 
192  m_AbsorberLogicalVolSet.insert(logic_absorber);
193 
194  G4LogicalVolume* logic_scint = new G4LogicalVolume(solid_scintillator,
195  material_scintillator,
196  "hHcal_scintillator_plate_logic",
197  0, 0, 0);
198  m_ScintiLogicalVolSet.insert(logic_scint);
199 
200  G4LogicalVolume* logic_wls = new G4LogicalVolume(solid_WLS_plate,
201  material_wls,
202  "hHcal_wls_plate_logic",
203  0, 0, 0);
204 
205  m_AbsorberLogicalVolSet.insert(logic_wls);
206  G4LogicalVolume* logic_support = new G4LogicalVolume(solid_support_plate,
207  material_support,
208  "hHcal_support_plate_logic",
209  0, 0, 0);
210 
211  m_AbsorberLogicalVolSet.insert(logic_support);
212  m_DisplayAction->AddVolume(logic_absorber, "Absorber");
213  m_DisplayAction->AddVolume(logic_scint, "Scintillator");
214  m_DisplayAction->AddVolume(logic_wls, "WLSplate");
215  m_DisplayAction->AddVolume(logic_support, "SupportPlate");
216 
217  /* place physical volumes for absorber and scintillator plates */
218  G4double xpos_i = -WlsDw / 2.0;
219  G4double ypos_i = -SupportDw / 2.0;
220  G4double zpos_i = (-1 * TowerDz / 2.0) + thickness_absorber / 2.0;
221 
222  string name_absorber = m_TowerLogicNamePrefix + "_single_plate_absorber";
223 
224  string name_scintillator = m_TowerLogicNamePrefix + "_single_plate_scintillator";
225 
226  string name_wls = m_TowerLogicNamePrefix + "_single_plate_wls";
227 
228  string name_support = m_TowerLogicNamePrefix + "_single_plate_support";
229 
230  for (int i = 1; i <= nlayers; i++)
231  {
232  new G4PVPlacement(0, G4ThreeVector(xpos_i, ypos_i, zpos_i),
233  logic_absorber,
234  name_absorber,
235  single_tower_logic,
236  0, 0, OverlapCheck());
237 
238  zpos_i += (thickness_absorber / 2. + thickness_scintillator / 2.);
239 
240  new G4PVPlacement(0, G4ThreeVector(xpos_i, ypos_i, zpos_i),
241  logic_scint,
242  name_scintillator,
243  single_tower_logic,
244  0, 0, OverlapCheck());
245 
246  zpos_i += (thickness_absorber / 2. + thickness_scintillator / 2.);
247  }
248  new G4PVPlacement(0, G4ThreeVector(0, (TowerDy / 2) - SupportDw / 2, 0),
249  logic_support,
250  name_support,
251  single_tower_logic,
252  0, 0, OverlapCheck());
253 
254  new G4PVPlacement(0, G4ThreeVector((TowerDx / 2) - WlsDw / 2, -SupportDw / 2, 0),
255  logic_wls,
256  name_wls,
257  single_tower_logic,
258  0, 0, OverlapCheck());
259  m_DisplayAction->AddVolume(single_tower_logic, "SingleTower");
260 
261  if (Verbosity() > 0)
262  {
263  std::cout << "PHG4BackwardHcalDetector: Building logical volume for single tower done." << std::endl;
264  }
265 
266  return single_tower_logic;
267 }
268 
270 {
271  /* Loop over all tower positions in vector and place tower */
272  for (std::map<std::string, towerposition>::iterator iterator = m_TowerPostionMap.begin(); iterator != m_TowerPostionMap.end(); ++iterator)
273  {
274  if (Verbosity() > 0)
275  {
276  std::cout << "PHG4BackwardHcalDetector: Place tower " << iterator->first
277  << " idx_j = " << iterator->second.idx_j << ", idx_k = " << iterator->second.idx_k
278  << " at x = " << iterator->second.x << " , y = " << iterator->second.y << " , z = " << iterator->second.z << std::endl;
279  }
280 
281  int copyno = (iterator->second.idx_j << 16) + iterator->second.idx_k;
282  new G4PVPlacement(0, G4ThreeVector(iterator->second.x, iterator->second.y, iterator->second.z),
283  singletower,
284  iterator->first,
285  hcalenvelope,
286  0, copyno, OverlapCheck());
287  }
288 
289  return 0;
290 }
291 
293 {
294  /* Open the datafile, if it won't open return an error */
295  ifstream istream_mapping;
296  istream_mapping.open(m_Params->get_string_param("mapping_file"));
297  if (!istream_mapping.is_open())
298  {
299  std::cout << "ERROR in PHG4BackwardHcalDetector: Failed to open mapping file " << m_Params->get_string_param("mapping_file") << std::endl;
300  gSystem->Exit(1);
301  }
302 
303  /* loop over lines in file */
304  string line_mapping;
305  while (getline(istream_mapping, line_mapping))
306  {
307  /* Skip lines starting with / including a '#' */
308  if (line_mapping.find("#") != string::npos)
309  {
310  if (Verbosity() > 0)
311  {
312  std::cout << "PHG4BackwardHcalDetector: SKIPPING line in mapping file: " << line_mapping << std::endl;
313  }
314  continue;
315  }
316 
317  istringstream iss(line_mapping);
318 
319  /* If line starts with keyword Tower, add to tower positions */
320  if (line_mapping.find("Tower ") != string::npos)
321  {
322  unsigned idx_j, idx_k, idx_l;
323  G4double pos_x, pos_y, pos_z;
324  G4double size_x, size_y, size_z;
325  G4double rot_x, rot_y, rot_z;
326  G4double dummy;
327  string dummys;
328 
329  /* read string- break if error */
330  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))
331  {
332  cout << "ERROR in PHG4BackwardHcalDetector: Failed to read line in mapping file " << m_Params->get_string_param("mapping_file") << std::endl;
333  gSystem->Exit(1);
334  }
335 
336  /* Construct unique name for tower */
337  /* Mapping file uses cm, this class uses mm for length */
338  ostringstream towername;
339  towername.str("");
340  towername << m_TowerLogicNamePrefix << "_j_" << idx_j << "_k_" << idx_k;
341 
342  /* Add Geant4 units */
343  pos_x = pos_x * cm;
344  pos_y = pos_y * cm;
345  pos_z = pos_z * cm;
346 
347  /* insert tower into tower map */
348  towerposition tower_new;
349  tower_new.x = pos_x;
350  tower_new.y = pos_y;
351  tower_new.z = pos_z;
352  tower_new.idx_j = idx_j;
353  tower_new.idx_k = idx_k;
354  m_TowerPostionMap.insert(make_pair(towername.str(), tower_new));
355  }
356  else
357  {
358  /* If this line is not a comment and not a tower, save parameter as string / value. */
359  string parname;
360  G4double parval;
361 
362  /* read string- break if error */
363  if (!(iss >> parname >> parval))
364  {
365  cout << "ERROR in PHG4BackwardHcalDetector: Failed to read line in mapping file " << m_Params->get_string_param("mapping_file") << std::endl;
366  gSystem->Exit(1);
367  }
368 
369  m_GlobalParameterMap.insert(make_pair(parname, parval));
370  }
371  }
372 
373  /* Update member variables for global parameters based on parsed parameter file */
374  std::map<string, G4double>::iterator parit;
375 
376  parit = m_GlobalParameterMap.find("Gtower_dx");
377  if (parit != m_GlobalParameterMap.end())
378  {
379  m_Params->set_double_param("tower_dx", parit->second); // in cm
380  }
381 
382  parit = m_GlobalParameterMap.find("Gtower_dy");
383  if (parit != m_GlobalParameterMap.end())
384  {
385  m_Params->set_double_param("tower_dy", parit->second); // in cm
386  }
387 
388  parit = m_GlobalParameterMap.find("Gtower_dz");
389  if (parit != m_GlobalParameterMap.end())
390  {
391  m_Params->set_double_param("tower_dz", parit->second); // in cm
392  }
393 
394  parit = m_GlobalParameterMap.find("Gr1_inner");
395  if (parit != m_GlobalParameterMap.end())
396  {
397  m_Params->set_double_param("rMin1", parit->second);
398  }
399 
400  parit = m_GlobalParameterMap.find("Gr1_outer");
401  if (parit != m_GlobalParameterMap.end())
402  {
403  m_Params->set_double_param("rMax1", parit->second);
404  }
405 
406  parit = m_GlobalParameterMap.find("Gr2_inner");
407  if (parit != m_GlobalParameterMap.end())
408  {
409  m_Params->set_double_param("rMin2", parit->second);
410  }
411 
412  parit = m_GlobalParameterMap.find("Gr2_outer");
413  if (parit != m_GlobalParameterMap.end())
414  {
415  m_Params->set_double_param("rMax2", parit->second);
416  }
417 
418  parit = m_GlobalParameterMap.find("Gdz");
419  if (parit != m_GlobalParameterMap.end())
420  {
421  m_Params->set_double_param("dZ", parit->second);
422  }
423 
424  parit = m_GlobalParameterMap.find("Gx0");
425  if (parit != m_GlobalParameterMap.end())
426  {
427  m_Params->set_double_param("place_x", parit->second);
428  }
429 
430  parit = m_GlobalParameterMap.find("Gy0");
431  if (parit != m_GlobalParameterMap.end())
432  {
433  m_Params->set_double_param("place_y", parit->second);
434  }
435 
436  parit = m_GlobalParameterMap.find("Gz0");
437  if (parit != m_GlobalParameterMap.end())
438  {
439  m_Params->set_double_param("place_z", parit->second);
440  }
441 
442  parit = m_GlobalParameterMap.find("Grot_x");
443  if (parit != m_GlobalParameterMap.end())
444  {
445  m_Params->set_double_param("rot_x", parit->second * rad / deg);
446  }
447 
448  parit = m_GlobalParameterMap.find("Grot_y");
449  if (parit != m_GlobalParameterMap.end())
450  {
451  m_Params->set_double_param("rot_y", parit->second * rad / deg);
452  }
453 
454  parit = m_GlobalParameterMap.find("Grot_z");
455  if (parit != m_GlobalParameterMap.end())
456  {
457  m_Params->set_double_param("rot_z", parit->second * rad / deg);
458  }
459 
460  return 0;
461 }